I need your help!

I want your feedback to make the book better for you and other readers. If you find typos, errors, or places where the text may be improved, please let me know. The best ways to provide feedback are by GitHub or hypothes.is annotations.

You can leave a comment at the bottom of the page/chapter, or open an issue or submit a pull request on GitHub: https://github.com/isaactpetersen/Fantasy-Football-Analytics-Textbook

Alternatively, you can leave an annotation using hypothes.is. To add an annotation, select some text and then click the symbol on the pop-up menu. To see the annotations of others, click the symbol in the upper right-hand corner of the page.

19  Machine Learning

19.1 Getting Started

19.1.1 Load Packages

library("petersenlab")
library("parallel")
library("future")
library("missRanger")
library("powerjoin")
library("tidymodels")
library("LongituRF")
library("gpboost")
library("tidyverse")

19.1.2 Load Data

# Downloaded Data - Processed
load(file = "./data/nfl_players.RData")
load(file = "./data/nfl_teams.RData")
load(file = "./data/nfl_rosters.RData")
load(file = "./data/nfl_rosters_weekly.RData")
load(file = "./data/nfl_schedules.RData")
load(file = "./data/nfl_combine.RData")
load(file = "./data/nfl_draftPicks.RData")
load(file = "./data/nfl_depthCharts.RData")
load(file = "./data/nfl_pbp.RData")
load(file = "./data/nfl_4thdown.RData")
load(file = "./data/nfl_participation.RData")
#load(file = "./data/nfl_actualFantasyPoints_weekly.RData")
load(file = "./data/nfl_injuries.RData")
load(file = "./data/nfl_snapCounts.RData")
load(file = "./data/nfl_espnQBR_seasonal.RData")
load(file = "./data/nfl_espnQBR_weekly.RData")
load(file = "./data/nfl_nextGenStats_weekly.RData")
load(file = "./data/nfl_advancedStatsPFR_seasonal.RData")
load(file = "./data/nfl_advancedStatsPFR_weekly.RData")
load(file = "./data/nfl_playerContracts.RData")
load(file = "./data/nfl_ftnCharting.RData")
load(file = "./data/nfl_playerIDs.RData")
load(file = "./data/nfl_rankings_draft.RData")
load(file = "./data/nfl_rankings_weekly.RData")
load(file = "./data/nfl_expectedFantasyPoints_weekly.RData")
load(file = "./data/nfl_expectedFantasyPoints_pbp.RData")

# Calculated Data - Processed
load(file = "./data/nfl_actualStats_career.RData")
load(file = "./data/nfl_actualStats_seasonal.RData")
load(file = "./data/player_stats_weekly.RData")
load(file = "./data/player_stats_seasonal.RData")

19.1.3 Specify Options

options(scipen = 999) # prevent scientific notation

19.2 Overview of Machine Learning

Machine learning takes us away from focusing on causal inference. Machine learning does not care about which processes are causal—i.e., which processes influence the outcome. Instead, machine learning cares about prediction—it cares about a predictor variable to the extent that it increases predictive accuracy regardless of whether it is causally related to the outcome.

Machine learning can be useful for leveraging big data and lots of predictor variable to develop predictive models with greater accuracy. However, many machine learning techniques are black boxes—it is often unclear how or why certain predictions are made, which can make it difficult to interpret the model’s decisions and understand the underlying relationships between variables. Machine learning tends to be a data-driven, atheoretical technique. This can result in overfitting. Thus, when estimating machine learning models, it is common to keep a hold-out sample for use in cross-validation to evaluate the extent of shrinkage of model coefficients. The data that the model is trained on is known as the “training data”. The data that the model was not trained on but is then is independently tested on—i.e., the hold-out sample—is the “test data”. Shrinkage occurs when predictor variables explain some random error variance in the original model. When the model is applied to an independent sample (i.e., the test data), the predictive model will likely not perform quite as well, and the regressions coefficients will tend to get smaller (i.e., shrink).

If the test data were collected as part of the same processes as the original data and were merely held out for purposes of analysis, this is called internal cross-validation. If the test data were collected separately from the original data used to train the model, this is called external cross-validation.

Most machine learning methods were developed with cross-sectional data in mind. That is, they assume that each person has only one observation on the outcome variable. However, with longitudinal data, each person has multiple observations on the outcome variable.

When performing machine learning, various approaches may help address this:

• transform data from long to wide form, so that each person has only one row
• when designing the training and test sets, keep all measurements from the same person in the same data object (either the training or test set); do not have some measurements from a given person in the training set and other measurements from the same person in the test set
• use a machine learning approach that accounts for the clustered/nested nature of the data

19.3 Types of Machine Learning

There are many approaches to machine learning. This chapter discusses several key ones:

• supervised learning
  • continuous outcome (i.e., regression)
    • linear regression
    • lasso regression
    • ridge regression
    • elastic net regression
  • categorical outcome (i.e., classification)
    • logistic regression
    • support vector machine
    • random forest
    • extreme gradient boosting
• unsupervised learning
  • clustering
  • principal component analysis
• semi-supervised learning
• reinforcement learning
  • deep learning
• ensemble

Ensemble machine learning methods combine multiple machine learning approaches with the goal that combining multiple approaches might lead to more accurate predictions that any one method might be able to achieve on its own.

19.3.1 Supervised Learning

[DEFINE SUPERVISED LEARNING]

Unlike linear and logistic regression, various machine learning techniques can handle multicollinearity, including LASSO regression, ridge regression, and elastic net regression. Least absolute shrinkage and selection operator (LASSO) regression helps perform selection of which predictor variables to keep in the model by shrinking some coefficients to zero. Ridge regression shrinks the coefficients of predictor variables toward zero, but not to zero, so it does not perform selection of which predictor variables to retain; this allows it to allow nonzero coefficients for multiple correlated predictor variables in the context of multicollinearity. Elastic net involves a combination of LASSO and ridge regression; it performs selection of which predictor variables to keep by shrinking the coefficients of some predictor variables to zero, and it shrinks the coefficients of some predictor variables toward zero, to address multicollinearity.

Unless interactions or nonlinear terms are specified, linear, logistic, LASSO, ridge, and elastic net regresstion do not account for interactions among the predictor variables or for nonlinear associations between the predictor variables and the outcome variable. By contrast, random forests and extreme gradient boosting do account for interactions among the predictor variables and for nonlinear associations between the predictor variables and the outcome variable.

19.3.2 Unsupervised Learning

[DEFINE UNSUPERVISED LEARNING]

We describe cluster analysis in Chapter 21. We describe principal component analysis in Chapter 23.

19.3.3 Semi-supervised Learning

[DEFINE SEMI-SUPERVISED LEARNING]

19.3.4 Reinforcement Learning

[DEFINE REINFORCEMENT LEARNING]

19.4 Data Processing

19.4.1 Prepare Data for Merging

# Prepare data for merging
#-todo: calculate years_of_experience
## Use common name for the same (gsis_id) ID variable

#nfl_actualFantasyPoints_player_weekly <- nfl_actualFantasyPoints_player_weekly %>% 
#  rename(gsis_id = player_id)
#
#nfl_actualFantasyPoints_player_seasonal <- nfl_actualFantasyPoints_player_seasonal %>% 
#  rename(gsis_id = player_id)

player_stats_seasonal_offense <- player_stats_seasonal %>% 
  filter(position_group %in% c("QB","RB","WR","TE")) %>% 
  rename(gsis_id = player_id)

player_stats_weekly_offense <- player_stats_weekly %>% 
  filter(position_group %in% c("QB","RB","WR","TE")) %>% 
  rename(gsis_id = player_id)

nfl_expectedFantasyPoints_weekly <- nfl_expectedFantasyPoints_weekly %>% 
  rename(gsis_id = player_id)

## Rename other variables to ensure common names

## Ensure variables with the same name have the same type
nfl_players <- nfl_players %>% 
  mutate(
    birth_date = as.Date(birth_date),
    jersey_number = as.character(jersey_number),
    gsis_it_id = as.character(gsis_it_id),
    years_of_experience = as.integer(years_of_experience))

player_stats_seasonal_offense <- player_stats_seasonal_offense %>% 
  mutate(
    birth_date = as.Date(birth_date),
    jersey_number = as.character(jersey_number),
    gsis_it_id = as.character(gsis_it_id))

nfl_rosters <- nfl_rosters %>% 
  mutate(
    draft_number = as.integer(draft_number))

nfl_rosters_weekly <- nfl_rosters_weekly %>% 
  mutate(
    draft_number = as.integer(draft_number))

nfl_depthCharts <- nfl_depthCharts %>% 
  mutate(
    season = as.integer(season))

nfl_expectedFantasyPoints_weekly <- nfl_expectedFantasyPoints_weekly %>% 
  mutate(
    season = as.integer(season),
    receptions = as.integer(receptions)) %>% 
  distinct(gsis_id, season, week, .keep_all = TRUE) # drop duplicated rows

## Rename variables
nfl_draftPicks <- nfl_draftPicks %>%
  rename(
    games_career = games,
    pass_completions_career = pass_completions,
    pass_attempts_career = pass_attempts,
    pass_yards_career = pass_yards,
    pass_tds_career = pass_tds,
    pass_ints_career = pass_ints,
    rush_atts_career = rush_atts,
    rush_yards_career = rush_yards,
    rush_tds_career = rush_tds,
    receptions_career = receptions,
    rec_yards_career = rec_yards,
    rec_tds_career = rec_tds,
    def_solo_tackles_career = def_solo_tackles,
    def_ints_career = def_ints,
    def_sacks_career = def_sacks
  )

## Subset variables
nfl_expectedFantasyPoints_weekly <- nfl_expectedFantasyPoints_weekly %>% 
  select(gsis_id:position, contains("_exp"), contains("_diff"), contains("_team")) #drop "raw stats" variables (e.g., rec_yards_gained) so they don't get coalesced with actual stats

# Check duplicate ids
player_stats_seasonal_offense %>% 
  group_by(gsis_id, season) %>% 
  filter(n() > 1) %>% 
  head()

nfl_advancedStatsPFR_seasonal %>% 
  group_by(gsis_id, season) %>% 
  filter(n() > 1, !is.na(gsis_id)) %>% 
  select(gsis_id, pfr_id, season, team, everything()) %>% 
  head()

Identify objects with shared variable names:

dplyr::intersect(
  names(nfl_players),
  names(nfl_draftPicks))

[1] "gsis_id"  "position"

length(na.omit(nfl_players$position)) # use by default (more cases)

[1] 21360

length(na.omit(nfl_draftPicks$position))

[1] 2855

dplyr::intersect(
  names(player_stats_seasonal_offense),
  names(nfl_advancedStatsPFR_seasonal))

[1] "gsis_id" "season"  "team"    "age"    

length(na.omit(player_stats_seasonal_offense$season)) # use by default (more cases)

[1] 14859

length(na.omit(nfl_advancedStatsPFR_seasonal$season))

[1] 10395

length(na.omit(player_stats_seasonal_offense$team)) # use by default (more cases)

[1] 14858

length(na.omit(nfl_advancedStatsPFR_seasonal$team))

[1] 10395

length(na.omit(player_stats_seasonal_offense$age)) # use by default (more cases)

[1] 14859

length(na.omit(nfl_advancedStatsPFR_seasonal$age))

[1] 10325

dplyr::intersect(
  names(nfl_rosters_weekly),
  names(nfl_expectedFantasyPoints_weekly))

[1] "gsis_id"   "season"    "week"      "position"  "full_name"

length(na.omit(nfl_rosters_weekly$season)) # use by default (more cases)

[1] 845134

length(na.omit(nfl_expectedFantasyPoints_weekly$season))

[1] 100272

length(na.omit(nfl_rosters_weekly$week)) # use by default (more cases)

[1] 841942

length(na.omit(nfl_expectedFantasyPoints_weekly$week))

[1] 100272

length(na.omit(nfl_rosters_weekly$position)) # use by default (more cases)

[1] 845101

length(na.omit(nfl_expectedFantasyPoints_weekly$position))

[1] 97815

length(na.omit(nfl_rosters_weekly$full_name)) # use by default (more cases)

[1] 845118

length(na.omit(nfl_expectedFantasyPoints_weekly$full_name))

[1] 97815

19.4.2 Merge Data

To merge data, we use the powerjoin package (Fabri, 2022):

# Create lists of objects to merge, depending on data structure: id; or id-season; or id-season-week
#-todo: remove redundant variables
playerListToMerge <- list(
  nfl_players %>% filter(!is.na(gsis_id)),
  nfl_draftPicks %>% filter(!is.na(gsis_id)) %>% select(-season)
)

playerSeasonListToMerge <- list(
  player_stats_seasonal_offense %>% filter(!is.na(gsis_id), !is.na(season)),
  nfl_advancedStatsPFR_seasonal %>% filter(!is.na(gsis_id), !is.na(season))
)

playerSeasonWeekListToMerge <- list(
  nfl_rosters_weekly %>% filter(!is.na(gsis_id), !is.na(season), !is.na(week)),
  #nfl_actualStats_offense_weekly,
  nfl_expectedFantasyPoints_weekly %>% filter(!is.na(gsis_id), !is.na(season), !is.na(week))
  #nfl_advancedStatsPFR_weekly,
)

playerSeasonWeekPositionListToMerge <- list(
  nfl_depthCharts %>% filter(!is.na(gsis_id), !is.na(season), !is.na(week))
)

# Merge data
playerMerged <- playerListToMerge %>% 
  reduce(
    powerjoin::power_full_join,
    by = c("gsis_id"),
    conflict = coalesce_xy) # where the objects have the same variable name (e.g., position), keep the values from object 1, unless it's NA, in which case use the relevant value from object 2

playerSeasonMerged <- playerSeasonListToMerge %>% 
  reduce(
    powerjoin::power_full_join,
    by = c("gsis_id","season"),
    conflict = coalesce_xy) # where the objects have the same variable name (e.g., team), keep the values from object 1, unless it's NA, in which case use the relevant value from object 2

playerSeasonWeekMerged <- playerSeasonWeekListToMerge %>% 
  reduce(
    powerjoin::power_full_join,
    by = c("gsis_id","season","week"),
    conflict = coalesce_xy) # where the objects have the same variable name (e.g., position), keep the values from object 1, unless it's NA, in which case use the relevant value from object 2

Identify objects with shared variable names:

dplyr::intersect(
  names(playerSeasonMerged),
  names(playerMerged))

 [1] "gsis_id"                  "position"                
 [3] "position_group"           "first_name"              
 [5] "last_name"                "esb_id"                  
 [7] "display_name"             "rookie_year"             
 [9] "college_conference"       "current_team_id"         
[11] "draft_club"               "draft_number"            
[13] "draftround"               "entry_year"              
[15] "football_name"            "gsis_it_id"              
[17] "headshot"                 "jersey_number"           
[19] "short_name"               "smart_id"                
[21] "status"                   "status_description_abbr" 
[23] "status_short_description" "uniform_number"          
[25] "height"                   "weight"                  
[27] "college_name"             "birth_date"              
[29] "suffix"                   "years_of_experience"     
[31] "pfr_player_name"          "team"                    
[33] "age"                     

seasonalData <- powerjoin::power_full_join(
  playerSeasonMerged,
  playerMerged %>% select(-age, -years_of_experience, -team, -team_abbr, -team_seq, -current_team_id), # drop variables from id objects that change from year to year (and thus are not necessarily accurate for a given season)
  by = "gsis_id",
  conflict = coalesce_xy # where the objects have the same variable name (e.g., position), keep the values from object 1, unless it's NA, in which case use the relevant value from object 2
) %>% 
  filter(!is.na(season)) %>% 
  select(gsis_id, season, player_display_name, position, team, games, everything())

dplyr::intersect(
  names(playerSeasonWeekMerged),
  names(seasonalData))

 [1] "gsis_id"                 "season"                 
 [3] "week"                    "team"                   
 [5] "jersey_number"           "status"                 
 [7] "first_name"              "last_name"              
 [9] "birth_date"              "height"                 
[11] "weight"                  "college"                
[13] "pfr_id"                  "headshot_url"           
[15] "status_description_abbr" "football_name"          
[17] "esb_id"                  "gsis_it_id"             
[19] "smart_id"                "entry_year"             
[21] "rookie_year"             "draft_club"             
[23] "draft_number"            "position"               

seasonalAndWeeklyData <- powerjoin::power_full_join(
  playerSeasonWeekMerged,
  seasonalData,
  by = c("gsis_id","season"),
  conflict = coalesce_xy # where the objects have the same variable name (e.g., position), keep the values from object 1, unless it's NA, in which case use the relevant value from object 2
) %>% 
  filter(!is.na(week)) %>% 
  select(gsis_id, season, week, full_name, position, team, everything())

# Duplicate cases
seasonalData %>% 
  group_by(gsis_id, season) %>% 
  filter(n() > 1) %>% 
  head()

seasonalAndWeeklyData %>% 
  group_by(gsis_id, season, week) %>% 
  filter(n() > 1) %>% 
  head()

19.4.3 Additional Processing

# Convert character and logical variables to factors
seasonalData <- seasonalData %>% 
  mutate(
    across(
      where(is.character),
      as.factor
    ),
    across(
      where(is.logical),
      as.factor
    )
  )

19.4.4 Fill in Missing Data for Static Variables

seasonalData <- seasonalData %>% 
  arrange(gsis_id, season) %>% 
  group_by(gsis_id) %>% 
  fill(
    player_name, player_display_name, pos, position, position_group,
    .direction = "downup") %>% 
  ungroup()

19.4.5 Create New Data Object for Merging with Later Predictions

newData_seasonal <- seasonalData %>% 
  filter(season == max(season, na.rm = TRUE))

19.4.6 Lag Fantasy Points

seasonalData_lag <- seasonalData %>% 
  arrange(gsis_id, season) %>% 
  group_by(gsis_id) %>% 
  mutate(
    fantasyPoints_lag = lead(fantasyPoints)
  ) %>% 
  ungroup()

seasonalData_lag %>% 
  select(gsis_id, player_display_name, season, fantasyPoints, fantasyPoints_lag) # verify that lagging worked as expected

19.4.7 Subset to Predictor Variables and Outcome Variable

seasonalData_lag %>% select_if(~class(.) == "Date")

seasonalData_lag %>% select_if(is.character)

seasonalData_lag %>% select_if(is.factor)

seasonalData_lag %>% select_if(is.logical)

dropVars <- c(
  "birth_date", "loaded", "full_name", "player_name", "player_display_name", "display_name", "suffix", "headshot_url", "player", "pos",
  "espn_id", "sportradar_id", "yahoo_id", "rotowire_id", "pff_id", "fantasy_data_id", "sleeper_id", "pfr_id",
  "pfr_player_id", "cfb_player_id", "pfr_player_name", "esb_id", "gsis_it_id", "smart_id",
  "college", "college_name", "team_abbr", "current_team_id", "college_conference", "draft_club", "status_description_abbr",
  "status_short_description", "short_name", "headshot", "uniform_number", "jersey_number", "first_name", "last_name",
  "football_name", "team")

seasonalData_lag_subset <- seasonalData_lag %>% 
  dplyr::select(-any_of(dropVars))

19.4.8 Separate by Position

seasonalData_lag_subsetQB <- seasonalData_lag_subset %>% 
  filter(position == "QB") %>% 
  select(
    gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic,
    height, weight, rookie_year, draft_number,
    fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag,
    completions:rushing_2pt_conversions, special_teams_tds, contains(".pass"), contains(".rush"))

seasonalData_lag_subsetRB <- seasonalData_lag_subset %>% 
  filter(position == "RB") %>% 
  select(
    gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic,
    height, weight, rookie_year, draft_number,
    fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag,
    carries:special_teams_tds, contains(".rush"), contains(".rec"))

seasonalData_lag_subsetWR <- seasonalData_lag_subset %>% 
  filter(position == "WR") %>% 
  select(
    gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic,
    height, weight, rookie_year, draft_number,
    fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag,
    carries:special_teams_tds, contains(".rush"), contains(".rec"))

seasonalData_lag_subsetTE <- seasonalData_lag_subset %>% 
  filter(position == "TE") %>% 
  select(
    gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic,
    height, weight, rookie_year, draft_number,
    fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag,
    carries:special_teams_tds, contains(".rush"), contains(".rec"))

19.4.9 Split into Test and Training Data

seasonalData_lag_qb_all <- seasonalData_lag_subsetQB
seasonalData_lag_rb_all <- seasonalData_lag_subsetRB
seasonalData_lag_wr_all <- seasonalData_lag_subsetWR
seasonalData_lag_te_all <- seasonalData_lag_subsetTE

set.seed(52242) # for reproducibility (to keep the same train/holdout players)

activeQBs <- unique(seasonalData_lag_qb_all$gsis_id[which(seasonalData_lag_qb_all$season == max(seasonalData_lag_qb_all$season, na.rm = TRUE))])
retiredQBs <- unique(seasonalData_lag_qb_all$gsis_id[which(seasonalData_lag_qb_all$gsis_id %ni% activeQBs)])
numQBs <- length(unique(seasonalData_lag_qb_all$gsis_id))
qbHoldoutIDs <- sample(retiredQBs, size = ceiling(.2 * numQBs)) # holdout 20% of players

activeRBs <- unique(seasonalData_lag_rb_all$gsis_id[which(seasonalData_lag_rb_all$season == max(seasonalData_lag_rb_all$season, na.rm = TRUE))])
retiredRBs <- unique(seasonalData_lag_rb_all$gsis_id[which(seasonalData_lag_rb_all$gsis_id %ni% activeRBs)])
numRBs <- length(unique(seasonalData_lag_rb_all$gsis_id))
rbHoldoutIDs <- sample(retiredRBs, size = ceiling(.2 * numRBs)) # holdout 20% of players

set.seed(52242) # for reproducibility (to keep the same train/holdout players); added here to prevent a downstream error with predict.missRanger() due to missingness; this suggests that an error can arise from including a player in the holdout sample who has missingness in particular variables; would be good to identify which player(s) in the holdout sample evoke that error to identify the kinds of missingness that yield the error

activeWRs <- unique(seasonalData_lag_wr_all$gsis_id[which(seasonalData_lag_wr_all$season == max(seasonalData_lag_wr_all$season, na.rm = TRUE))])
retiredWRs <- unique(seasonalData_lag_wr_all$gsis_id[which(seasonalData_lag_wr_all$gsis_id %ni% activeWRs)])
numWRs <- length(unique(seasonalData_lag_wr_all$gsis_id))
wrHoldoutIDs <- sample(retiredWRs, size = ceiling(.2 * numWRs)) # holdout 20% of players

activeTEs <- unique(seasonalData_lag_te_all$gsis_id[which(seasonalData_lag_te_all$season == max(seasonalData_lag_te_all$season, na.rm = TRUE))])
retiredTEs <- unique(seasonalData_lag_te_all$gsis_id[which(seasonalData_lag_te_all$gsis_id %ni% activeTEs)])
numTEs <- length(unique(seasonalData_lag_te_all$gsis_id))
teHoldoutIDs <- sample(retiredTEs, size = ceiling(.2 * numTEs)) # holdout 20% of players
  
seasonalData_lag_qb_train <- seasonalData_lag_qb_all %>% 
  filter(gsis_id %ni% qbHoldoutIDs)
seasonalData_lag_qb_test <- seasonalData_lag_qb_all %>% 
  filter(gsis_id %in% qbHoldoutIDs)

seasonalData_lag_rb_train <- seasonalData_lag_rb_all %>% 
  filter(gsis_id %ni% rbHoldoutIDs)
seasonalData_lag_rb_test <- seasonalData_lag_rb_all %>% 
  filter(gsis_id %in% rbHoldoutIDs)

seasonalData_lag_wr_train <- seasonalData_lag_wr_all %>% 
  filter(gsis_id %ni% wrHoldoutIDs)
seasonalData_lag_wr_test <- seasonalData_lag_wr_all %>% 
  filter(gsis_id %in% wrHoldoutIDs)

seasonalData_lag_te_train <- seasonalData_lag_te_all %>% 
  filter(gsis_id %ni% teHoldoutIDs)
seasonalData_lag_te_test <- seasonalData_lag_te_all %>% 
  filter(gsis_id %in% teHoldoutIDs)

19.4.10 Impute the Missing Data

Here is a vignette demonstrating how to impute missing data using missForest(): https://rpubs.com/lmorgan95/MissForest (archived at: https://perma.cc/6GB4-2E22). Below, we impute the training data (and all data) separately by position. We then use the imputed training data to make out-of-sample predictions to fill in the missing data for the testing data. We do not want to impute the training and testing data together so that we can keep them separate for the purposes of cross-validation. However, we impute all data (training and test data together) for purposes of making out-of-sample predictions from the machine learning models to predict players’ performance next season (when actuals are not yet available for evaluating their accuracy). To impute data, we use the missRanger package (Mayer, 2024).

Note 19.1: Impute missing data for machine learning

Note: the following code takes a while to run.

# QBs
seasonalData_lag_qb_all_imp <- missRanger::missRanger(
  seasonalData_lag_qb_all,
  pmm.k = 5,
  verbose = 2,
  seed = 52242,
  keep_forests = TRUE)

Variables to impute:        fantasy_points, fantasy_points_ppr, special_teams_tds, passing_epa, pacr, rushing_epa, fantasyPoints_lag, passing_cpoe, rookie_year, draft_number, gs, pass_attempts.pass, throwaways.pass, spikes.pass, drops.pass, bad_throws.pass, times_blitzed.pass, times_hurried.pass, times_hit.pass, times_pressured.pass, batted_balls.pass, on_tgt_throws.pass, rpo_plays.pass, rpo_yards.pass, rpo_pass_att.pass, rpo_pass_yards.pass, rpo_rush_att.pass, rpo_rush_yards.pass, pa_pass_att.pass, pa_pass_yards.pass, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, drop_pct.pass, bad_throw_pct.pass, on_tgt_pct.pass, pressure_pct.pass, ybc_att.rush, yac_att.rush, pocket_time.pass
Variables used to impute:   gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic, height, weight, rookie_year, draft_number, fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag, completions, attempts, passing_yards, passing_tds, passing_interceptions, sacks_suffered, sack_yards_lost, sack_fumbles, sack_fumbles_lost, passing_air_yards, passing_yards_after_catch, passing_first_downs, passing_epa, passing_cpoe, passing_2pt_conversions, pacr, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_epa, rushing_2pt_conversions, special_teams_tds, pocket_time.pass, pass_attempts.pass, throwaways.pass, spikes.pass, drops.pass, bad_throws.pass, times_blitzed.pass, times_hurried.pass, times_hit.pass, times_pressured.pass, batted_balls.pass, on_tgt_throws.pass, rpo_plays.pass, rpo_yards.pass, rpo_pass_att.pass, rpo_pass_yards.pass, rpo_rush_att.pass, rpo_rush_yards.pass, pa_pass_att.pass, pa_pass_yards.pass, drop_pct.pass, bad_throw_pct.pass, on_tgt_pct.pass, pressure_pct.pass, ybc_att.rush, yac_att.rush, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush

    fntsy_  fnts__  spcl__  pssng_p pacr    rshng_  fntsP_  pssng_c rok_yr  drft_n  gs  pss_t.  thrww.  spks.p  drps.p  bd_th.  tms_b.  tms_hr. tms_ht. tms_p.  bttd_.  on_tgt_t.   rp_pl.  rp_yr.  rp_pss_t.   rp_pss_y.   rp_rsh_t.   rp_rsh_y.   p_pss_t.    p_pss_y.    att.rs  yds.rs  td.rsh  x1d.rs  ybc.rs  yc.rsh  brk_t.  att_b.  drp_p.  bd_t_.  on_tgt_p.   prss_.  ybc_t.  yc_tt.  pckt_.
iter 1: 0.0054  0.0024  0.7924  0.1919  0.7612  0.3628  0.4789  0.4133  0.0224  0.5216  0.0271  0.0134  0.3024  0.7659  0.1304  0.0541  0.0758  0.1759  0.1820  0.0370  0.3238  0.0291  0.2952  0.1812  0.0885  0.0867  0.2627  0.2563  0.1093  0.0902  0.0580  0.0645  0.1732  0.0524  0.0578  0.1795  0.3524  0.3428  0.7447  0.5158  0.0824  0.6803  0.3529  0.5758  0.8111  
iter 2: 0.0044  0.0048  0.8304  0.2002  0.7926  0.3736  0.4801  0.4289  0.0488  0.6139  0.0188  0.0090  0.2883  0.7481  0.0764  0.0385  0.0718  0.1231  0.1329  0.0337  0.2760  0.0113  0.0548  0.0814  0.0765  0.0990  0.1989  0.2841  0.0707  0.0952  0.0396  0.0386  0.1606  0.0492  0.0525  0.1220  0.2541  0.3556  0.7468  0.4937  0.0827  0.6610  0.3465  0.5796  0.8134  
iter 3: 0.0049  0.0046  0.8690  0.1986  0.7810  0.3641  0.4774  0.4360  0.0528  0.6123  0.0188  0.0088  0.2867  0.7538  0.0767  0.0393  0.0734  0.1261  0.1374  0.0343  0.2741  0.0119  0.0524  0.0816  0.0748  0.1008  0.2184  0.2811  0.0691  0.0926  0.0389  0.0413  0.1640  0.0511  0.0585  0.1255  0.2510  0.3609  0.7477  0.5108  0.0858  0.6426  0.3588  0.5734  0.8300  

seasonalData_lag_qb_all_imp

missRanger object. Extract imputed data via $data
- best iteration: 2 
- best average OOB imputation error: 0.2524825 

data_all_qb <- seasonalData_lag_qb_all_imp$data
data_all_qb$fantasyPointsMC_lag <- scale(data_all_qb$fantasyPoints_lag, scale = FALSE) # mean-centered
data_all_qb_matrix <- data_all_qb %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()
newData_qb <- data_all_qb %>% 
  filter(season == max(season, na.rm = TRUE)) %>% 
  select(-fantasyPoints_lag, -fantasyPointsMC_lag)
newData_qb_matrix <- data_all_qb_matrix[
  data_all_qb_matrix[, "season"] == max(data_all_qb_matrix[, "season"], na.rm = TRUE), # keep only rows with the most recent season
  , # all columns
  drop = FALSE]

dropCol_qb <- which(colnames(newData_qb_matrix) %in% c("fantasyPoints_lag","fantasyPointsMC_lag"))
newData_qb_matrix <- newData_qb_matrix[, -dropCol_qb, drop = FALSE]

seasonalData_lag_qb_train_imp <- missRanger::missRanger(
  seasonalData_lag_qb_train,
  pmm.k = 5,
  verbose = 2,
  seed = 52242,
  keep_forests = TRUE)

Variables to impute:        fantasy_points, fantasy_points_ppr, special_teams_tds, passing_epa, pacr, rushing_epa, fantasyPoints_lag, passing_cpoe, rookie_year, draft_number, gs, pass_attempts.pass, throwaways.pass, spikes.pass, drops.pass, bad_throws.pass, times_blitzed.pass, times_hurried.pass, times_hit.pass, times_pressured.pass, batted_balls.pass, on_tgt_throws.pass, rpo_plays.pass, rpo_yards.pass, rpo_pass_att.pass, rpo_pass_yards.pass, rpo_rush_att.pass, rpo_rush_yards.pass, pa_pass_att.pass, pa_pass_yards.pass, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, drop_pct.pass, bad_throw_pct.pass, on_tgt_pct.pass, pressure_pct.pass, ybc_att.rush, yac_att.rush, pocket_time.pass
Variables used to impute:   gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic, height, weight, rookie_year, draft_number, fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag, completions, attempts, passing_yards, passing_tds, passing_interceptions, sacks_suffered, sack_yards_lost, sack_fumbles, sack_fumbles_lost, passing_air_yards, passing_yards_after_catch, passing_first_downs, passing_epa, passing_cpoe, passing_2pt_conversions, pacr, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_epa, rushing_2pt_conversions, special_teams_tds, pocket_time.pass, pass_attempts.pass, throwaways.pass, spikes.pass, drops.pass, bad_throws.pass, times_blitzed.pass, times_hurried.pass, times_hit.pass, times_pressured.pass, batted_balls.pass, on_tgt_throws.pass, rpo_plays.pass, rpo_yards.pass, rpo_pass_att.pass, rpo_pass_yards.pass, rpo_rush_att.pass, rpo_rush_yards.pass, pa_pass_att.pass, pa_pass_yards.pass, drop_pct.pass, bad_throw_pct.pass, on_tgt_pct.pass, pressure_pct.pass, ybc_att.rush, yac_att.rush, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush

    fntsy_  fnts__  spcl__  pssng_p pacr    rshng_  fntsP_  pssng_c rok_yr  drft_n  gs  pss_t.  thrww.  spks.p  drps.p  bd_th.  tms_b.  tms_hr. tms_ht. tms_p.  bttd_.  on_tgt_t.   rp_pl.  rp_yr.  rp_pss_t.   rp_pss_y.   rp_rsh_t.   rp_rsh_y.   p_pss_t.    p_pss_y.    att.rs  yds.rs  td.rsh  x1d.rs  ybc.rs  yc.rsh  brk_t.  att_b.  drp_p.  bd_t_.  on_tgt_p.   prss_.  ybc_t.  yc_tt.  pckt_.
iter 1: 0.0061  0.0028  0.8162  0.1897  0.5083  0.3633  0.4726  0.4456  0.0242  0.4723  0.0283  0.0141  0.2939  0.7728  0.1343  0.0558  0.0744  0.1757  0.1818  0.0381  0.3288  0.0351  0.2921  0.1846  0.0860  0.0894  0.2737  0.2661  0.1127  0.0900  0.0586  0.0644  0.1800  0.0574  0.0639  0.1792  0.3570  0.3486  0.7646  0.5313  0.0868  0.7084  0.3533  0.5933  0.8466  
iter 2: 0.0052  0.0052  0.8304  0.1937  0.5621  0.3715  0.4614  0.4586  0.0505  0.5647  0.0192  0.0092  0.2953  0.7530  0.0800  0.0393  0.0725  0.1170  0.1355  0.0343  0.2771  0.0121  0.0555  0.0731  0.0713  0.0979  0.2073  0.2943  0.0698  0.0911  0.0416  0.0399  0.1683  0.0527  0.0577  0.1262  0.2474  0.3582  0.7719  0.5165  0.0900  0.6862  0.3642  0.5926  0.8400  
iter 3: 0.0053  0.0051  0.8261  0.2008  0.5551  0.3571  0.4727  0.4410  0.0551  0.5658  0.0188  0.0092  0.2859  0.7460  0.0807  0.0402  0.0739  0.1202  0.1393  0.0351  0.2808  0.0114  0.0595  0.0705  0.0775  0.1051  0.2163  0.2935  0.0718  0.0921  0.0426  0.0400  0.1719  0.0535  0.0534  0.1225  0.2498  0.3484  0.7502  0.5100  0.0884  0.6609  0.3672  0.5852  0.8440  
iter 4: 0.0054  0.0051  0.6928  0.1979  0.5598  0.3732  0.4771  0.4349  0.0506  0.5691  0.0189  0.0085  0.2891  0.7456  0.0785  0.0395  0.0737  0.1210  0.1353  0.0335  0.2836  0.0117  0.0566  0.0778  0.0743  0.1055  0.2131  0.2964  0.0697  0.0912  0.0396  0.0395  0.1611  0.0531  0.0597  0.1258  0.2600  0.3560  0.8062  0.5032  0.0973  0.6739  0.3698  0.5875  0.8485  
iter 5: 0.0052  0.0055  0.8355  0.1965  0.5664  0.3710  0.4743  0.4604  0.0520  0.5598  0.0193  0.0091  0.2852  0.7474  0.0800  0.0405  0.0722  0.1213  0.1366  0.0344  0.2788  0.0118  0.0555  0.0756  0.0746  0.0986  0.2190  0.2765  0.0695  0.0932  0.0390  0.0425  0.1650  0.0509  0.0576  0.1305  0.2556  0.3509  0.7738  0.5051  0.0969  0.6902  0.3640  0.6007  0.8326  

seasonalData_lag_qb_train_imp

missRanger object. Extract imputed data via $data
- best iteration: 4 
- best average OOB imputation error: 0.2482278 

data_train_qb <- seasonalData_lag_qb_train_imp$data
data_train_qb$fantasyPointsMC_lag <- scale(data_train_qb$fantasyPoints_lag, scale = FALSE) # mean-centered
data_train_qb_matrix <- data_train_qb %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()

seasonalData_lag_qb_test_imp <- predict(
  object = seasonalData_lag_qb_train_imp,
  newdata = seasonalData_lag_qb_test,
  seed = 52242)

data_test_qb <- seasonalData_lag_qb_test_imp
data_test_qb_matrix <- data_test_qb %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()

# RBs
seasonalData_lag_rb_all_imp <- missRanger::missRanger(
  seasonalData_lag_rb_all,
  pmm.k = 5,
  verbose = 2,
  seed = 52242,
  keep_forests = TRUE)

Variables to impute:        games, ageCentered20, ageCentered20Quadratic, fantasy_points, fantasy_points_ppr, fantasyPoints, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_2pt_conversions, receptions, targets, receiving_yards, receiving_tds, receiving_fumbles, receiving_fumbles_lost, receiving_air_yards, receiving_yards_after_catch, receiving_first_downs, receiving_2pt_conversions, special_teams_tds, years_of_experience, rushing_epa, air_yards_share, receiving_epa, racr, target_share, wopr, fantasyPoints_lag, rookie_year, draft_number, gs, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, ybc_att.rush, yac_att.rush, adot.rec, rat.rec, drop_percent.rec, rec_br.rec, ybc_r.rec, yac_r.rec
Variables used to impute:   gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic, height, weight, rookie_year, draft_number, fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_epa, rushing_2pt_conversions, receptions, targets, receiving_yards, receiving_tds, receiving_fumbles, receiving_fumbles_lost, receiving_air_yards, receiving_yards_after_catch, receiving_first_downs, receiving_epa, receiving_2pt_conversions, racr, target_share, air_yards_share, wopr, special_teams_tds, ybc_att.rush, yac_att.rush, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, ybc_r.rec, yac_r.rec, adot.rec, rat.rec, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, drop_percent.rec, rec_br.rec

    games   agCn20  agC20Q  fntsy_  fnts__  fntsyP  carris  rshng_y rshng_t rshng_f rshng_fm_   rshng_fr_   rsh_2_  rcptns  targts  rcvng_y rcvng_t rcvng_f rcvng_fm_   rcvng_r_    rcv___  rcvng_fr_   rcv_2_  spcl__  yrs_f_  rshng_p ar_yr_  rcvng_p racr    trgt_s  wopr    fntsP_  rok_yr  drft_n  gs  att.rs  yds.rs  td.rsh  x1d.rs  ybc.rs  yc.rsh  brk_tkl.rs  att_b.  tgt.rc  rec.rc  yds.rc  td.rec  x1d.rc  ybc.rc  yac.rc  brk_tkl.rc  drp.rc  int.rc  ybc_t.  yc_tt.  adt.rc  rat.rc  drp_p.  rc_br.  ybc_r.  yc_r.r
iter 1: 0.8865  0.0057  0.0031  0.4544  0.0178  0.0032  0.0745  0.0233  0.1462  0.4895  0.2594  0.0295  0.9849  0.0690  0.0666  0.0534  0.4327  0.8626  0.4824  0.6841  0.0322  0.0614  1.0171  0.8263  0.1817  0.4512  0.3321  0.3894  0.5211  0.4549  0.1817  0.5440  0.0197  0.5999  0.1700  0.0244  0.0222  0.0792  0.0297  0.0527  0.0520  0.2134  0.3431  0.0252  0.0180  0.0257  0.1634  0.0437  0.3108  0.0217  0.3925  0.4610  0.6941  0.4880  0.5402  0.2670  0.2026  0.3482  0.1596  0.2698  0.3637  
iter 2: 0.2755  0.0162  0.0207  0.0063  0.0037  0.0044  0.0161  0.0148  0.0912  0.2524  0.2898  0.0248  0.9832  0.0273  0.0444  0.0233  0.2065  0.4600  0.4891  0.1285  0.0332  0.0457  1.0175  0.8566  0.1824  0.4212  0.2329  0.3099  0.5605  0.2569  0.1742  0.5373  0.0424  0.6377  0.1653  0.0167  0.0123  0.0859  0.0302  0.0367  0.0349  0.1030  0.3689  0.0144  0.0159  0.0195  0.1403  0.0434  0.1549  0.0190  0.3840  0.1050  0.5453  0.4882  0.5616  0.2472  0.1953  0.1525  0.1687  0.2595  0.3619  
iter 3: 0.2744  0.0163  0.0231  0.0062  0.0038  0.0047  0.0152  0.0137  0.0980  0.2601  0.2906  0.0244  0.9800  0.0265  0.0347  0.0231  0.2101  0.4638  0.4954  0.1284  0.0283  0.0458  1.0114  0.8731  0.1818  0.4117  0.2278  0.3037  0.5699  0.2052  0.1800  0.5389  0.0400  0.6423  0.1624  0.0166  0.0124  0.0893  0.0306  0.0374  0.0356  0.1074  0.3628  0.0144  0.0163  0.0187  0.1390  0.0463  0.1583  0.0190  0.3882  0.1062  0.5642  0.4796  0.5570  0.2380  0.1935  0.1586  0.1588  0.2625  0.3648  
iter 4: 0.2776  0.0169  0.0220  0.0063  0.0038  0.0045  0.0151  0.0138  0.0979  0.2584  0.2846  0.0243  0.9782  0.0263  0.0281  0.0221  0.1968  0.4594  0.4817  0.1267  0.0290  0.0462  1.0104  0.8614  0.1854  0.4216  0.2333  0.3004  0.5467  0.1917  0.1815  0.5353  0.0443  0.6503  0.1657  0.0166  0.0121  0.0905  0.0313  0.0378  0.0357  0.1041  0.3437  0.0155  0.0159  0.0185  0.1405  0.0441  0.1613  0.0196  0.3816  0.1117  0.5682  0.5011  0.5585  0.2421  0.1975  0.1520  0.1770  0.2650  0.3647  
iter 5: 0.2752  0.0163  0.0226  0.0063  0.0038  0.0045  0.0158  0.0138  0.1015  0.2614  0.2857  0.0242  0.9740  0.0250  0.0303  0.0218  0.2004  0.4607  0.4810  0.1167  0.0285  0.0449  1.0077  0.8658  0.1835  0.4182  0.2170  0.2995  0.5690  0.2010  0.1794  0.5375  0.0385  0.6487  0.1652  0.0166  0.0124  0.0878  0.0306  0.0368  0.0353  0.1069  0.3539  0.0154  0.0159  0.0193  0.1409  0.0447  0.1598  0.0205  0.3873  0.1062  0.5583  0.4895  0.5501  0.2418  0.1979  0.1713  0.1726  0.2625  0.3596  
iter 6: 0.2760  0.0158  0.0223  0.0063  0.0037  0.0046  0.0150  0.0144  0.0982  0.2568  0.2816  0.0238  0.9810  0.0253  0.0273  0.0223  0.2141  0.4606  0.4881  0.1386  0.0300  0.0457  1.0174  0.8605  0.1821  0.4188  0.2263  0.2985  0.5497  0.1779  0.1536  0.5388  0.0389  0.6422  0.1668  0.0162  0.0119  0.0897  0.0305  0.0376  0.0356  0.1066  0.3529  0.0149  0.0159  0.0196  0.1446  0.0450  0.1585  0.0197  0.3857  0.1001  0.5607  0.4948  0.5478  0.2487  0.1945  0.1438  0.1543  0.2568  0.3612  
iter 7: 0.2748  0.0158  0.0212  0.0064  0.0039  0.0047  0.0149  0.0141  0.0986  0.2611  0.2877  0.0241  0.9755  0.0253  0.0310  0.0223  0.2163  0.4553  0.4885  0.1335  0.0293  0.0456  1.0096  0.8575  0.1821  0.4236  0.2203  0.2998  0.5510  0.2107  0.1797  0.5354  0.0416  0.6395  0.1646  0.0166  0.0117  0.0895  0.0310  0.0371  0.0361  0.1073  0.3547  0.0154  0.0156  0.0193  0.1410  0.0449  0.1628  0.0201  0.3892  0.1076  0.5609  0.4946  0.5643  0.2392  0.1899  0.1540  0.1423  0.2667  0.3603  

seasonalData_lag_rb_all_imp

missRanger object. Extract imputed data via $data
- best iteration: 6 
- best average OOB imputation error: 0.2175486 

data_all_rb <- seasonalData_lag_rb_all_imp$data
data_all_rb$fantasyPointsMC_lag <- scale(data_all_rb$fantasyPoints_lag, scale = FALSE) # mean-centered
data_all_rb_matrix <- data_all_rb %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()
newData_rb <- data_all_rb %>% 
  filter(season == max(season, na.rm = TRUE)) %>% 
  select(-fantasyPoints_lag, -fantasyPointsMC_lag)
newData_rb_matrix <- data_all_rb_matrix[
  data_all_rb_matrix[, "season"] == max(data_all_rb_matrix[, "season"], na.rm = TRUE), # keep only rows with the most recent season
  , # all columns
  drop = FALSE]

dropCol_rb <- which(colnames(newData_rb_matrix) %in% c("fantasyPoints_lag","fantasyPointsMC_lag"))
newData_rb_matrix <- newData_rb_matrix[, -dropCol_rb, drop = FALSE]

seasonalData_lag_rb_train_imp <- missRanger::missRanger(
  seasonalData_lag_rb_train,
  pmm.k = 5,
  verbose = 2,
  seed = 52242,
  keep_forests = TRUE)

Variables to impute:        games, ageCentered20, ageCentered20Quadratic, fantasy_points, fantasy_points_ppr, fantasyPoints, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_2pt_conversions, receptions, targets, receiving_yards, receiving_tds, receiving_fumbles, receiving_fumbles_lost, receiving_air_yards, receiving_yards_after_catch, receiving_first_downs, receiving_2pt_conversions, special_teams_tds, years_of_experience, rushing_epa, air_yards_share, receiving_epa, racr, target_share, wopr, fantasyPoints_lag, rookie_year, draft_number, gs, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, ybc_att.rush, yac_att.rush, adot.rec, rat.rec, drop_percent.rec, rec_br.rec, ybc_r.rec, yac_r.rec
Variables used to impute:   gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic, height, weight, rookie_year, draft_number, fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_epa, rushing_2pt_conversions, receptions, targets, receiving_yards, receiving_tds, receiving_fumbles, receiving_fumbles_lost, receiving_air_yards, receiving_yards_after_catch, receiving_first_downs, receiving_epa, receiving_2pt_conversions, racr, target_share, air_yards_share, wopr, special_teams_tds, ybc_att.rush, yac_att.rush, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, ybc_r.rec, yac_r.rec, adot.rec, rat.rec, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, drop_percent.rec, rec_br.rec

    games   agCn20  agC20Q  fntsy_  fnts__  fntsyP  carris  rshng_y rshng_t rshng_f rshng_fm_   rshng_fr_   rsh_2_  rcptns  targts  rcvng_y rcvng_t rcvng_f rcvng_fm_   rcvng_r_    rcv___  rcvng_fr_   rcv_2_  spcl__  yrs_f_  rshng_p ar_yr_  rcvng_p racr    trgt_s  wopr    fntsP_  rok_yr  drft_n  gs  att.rs  yds.rs  td.rsh  x1d.rs  ybc.rs  yc.rsh  brk_tkl.rs  att_b.  tgt.rc  rec.rc  yds.rc  td.rec  x1d.rc  ybc.rc  yac.rc  brk_tkl.rc  drp.rc  int.rc  ybc_t.  yc_tt.  adt.rc  rat.rc  drp_p.  rc_br.  ybc_r.  yc_r.r
iter 1: 0.8759  0.0072  0.0036  0.4578  0.0178  0.0035  0.0736  0.0229  0.1524  0.4776  0.2679  0.0288  0.9965  0.0749  0.0744  0.0553  0.4578  0.8604  0.4998  0.6821  0.0360  0.0639  1.0042  0.8380  0.1806  0.4662  0.3419  0.3961  0.5595  0.4715  0.1968  0.5338  0.0246  0.5882  0.1726  0.0265  0.0235  0.0849  0.0311  0.0550  0.0521  0.2131  0.3689  0.0281  0.0197  0.0286  0.1742  0.0463  0.3114  0.0229  0.3942  0.4806  0.7239  0.5199  0.5631  0.2865  0.2195  0.3630  0.2052  0.2596  0.4091  
iter 2: 0.2745  0.0177  0.0266  0.0067  0.0041  0.0049  0.0169  0.0154  0.1017  0.2590  0.2956  0.0237  0.9814  0.0286  0.0522  0.0240  0.2187  0.4582  0.4919  0.1541  0.0362  0.0475  1.0075  0.8811  0.1822  0.4481  0.2377  0.3184  0.6116  0.2628  0.2007  0.5254  0.0473  0.6411  0.1653  0.0179  0.0132  0.0940  0.0325  0.0392  0.0375  0.1057  0.3678  0.0161  0.0169  0.0196  0.1510  0.0484  0.1521  0.0202  0.3941  0.1035  0.5567  0.5120  0.5666  0.2466  0.2087  0.1733  0.1699  0.2524  0.4066  
iter 3: 0.2766  0.0190  0.0273  0.0067  0.0041  0.0048  0.0159  0.0149  0.0971  0.2615  0.2952  0.0245  0.9668  0.0278  0.0409  0.0240  0.2180  0.4648  0.4931  0.1319  0.0350  0.0495  1.0128  0.8907  0.1820  0.4366  0.2459  0.3124  0.6236  0.2555  0.2114  0.5276  0.0438  0.6314  0.1658  0.0175  0.0122  0.0899  0.0319  0.0386  0.0386  0.1100  0.3783  0.0155  0.0165  0.0194  0.1477  0.0474  0.1499  0.0194  0.3929  0.1121  0.5761  0.5245  0.5651  0.2490  0.2103  0.1767  0.1817  0.2658  0.4106  

seasonalData_lag_rb_train_imp

missRanger object. Extract imputed data via $data
- best iteration: 2 
- best average OOB imputation error: 0.226086 

data_train_rb <- seasonalData_lag_rb_train_imp$data
data_train_rb$fantasyPointsMC_lag <- scale(data_train_rb$fantasyPoints_lag, scale = FALSE) # mean-centered
data_train_rb_matrix <- data_train_rb %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()

seasonalData_lag_rb_test_imp <- predict(
  object = seasonalData_lag_rb_train_imp,
  newdata = seasonalData_lag_rb_test,
  seed = 52242)

data_test_rb <- seasonalData_lag_rb_test_imp
data_test_rb_matrix <- data_test_rb %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()

# WRs
seasonalData_lag_wr_all_imp <- missRanger::missRanger(
  seasonalData_lag_wr_all,
  pmm.k = 5,
  verbose = 2,
  seed = 52242,
  keep_forests = TRUE)

Variables to impute:        fantasy_points, fantasy_points_ppr, special_teams_tds, years_of_experience, receiving_epa, racr, air_yards_share, target_share, wopr, fantasyPoints_lag, rookie_year, rushing_epa, draft_number, gs, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, adot.rec, rat.rec, drop_percent.rec, rec_br.rec, ybc_r.rec, yac_r.rec, ybc_att.rush, yac_att.rush
Variables used to impute:   gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic, height, weight, rookie_year, draft_number, fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_epa, rushing_2pt_conversions, receptions, targets, receiving_yards, receiving_tds, receiving_fumbles, receiving_fumbles_lost, receiving_air_yards, receiving_yards_after_catch, receiving_first_downs, receiving_epa, receiving_2pt_conversions, racr, target_share, air_yards_share, wopr, special_teams_tds, ybc_att.rush, yac_att.rush, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, ybc_r.rec, yac_r.rec, adot.rec, rat.rec, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, drop_percent.rec, rec_br.rec

    fntsy_  fnts__  spcl__  yrs_f_  rcvng_  racr    ar_yr_  trgt_s  wopr    fntsP_  rok_yr  rshng_  drft_n  gs  att.rs  yds.rs  td.rsh  x1d.rs  ybc.rs  yc.rsh  brk_tkl.rs  att_b.  tgt.rc  rec.rc  yds.rc  td.rec  x1d.rc  ybc.rc  yac.rc  brk_tkl.rc  drp.rc  int.rc  adt.rc  rat.rc  drp_p.  rc_br.  ybc_r.  yc_r.r  ybc_t.  yc_tt.
iter 1: 0.0061  0.0010  0.7104  0.1566  0.1040  0.8131  0.1013  0.1722  0.0402  0.4890  0.0150  0.3811  0.6654  0.1459  0.1184  0.0898  0.2353  0.1234  0.0966  0.2670  0.6383  0.3084  0.0198  0.0136  0.0151  0.0671  0.0135  0.0268  0.0442  0.4465  0.4320  0.4674  0.2961  0.1410  0.3819  0.1840  0.2251  0.3929  0.2568  0.4760  
iter 2: 0.0058  0.0019  0.7826  0.1601  0.0835  0.7518  0.0607  0.0930  0.0452  0.4939  0.0296  0.3301  0.6843  0.1440  0.0851  0.0600  0.2638  0.1161  0.0708  0.1804  0.3103  0.3223  0.0109  0.0108  0.0096  0.0719  0.0139  0.0200  0.0318  0.4476  0.0778  0.3692  0.2401  0.1448  0.1629  0.1601  0.2261  0.3793  0.2536  0.4775  
iter 3: 0.0061  0.0019  0.7857  0.1593  0.0829  0.7421  0.0580  0.0986  0.0481  0.4946  0.0318  0.3334  0.6890  0.1430  0.0823  0.0604  0.2595  0.1177  0.0728  0.1802  0.3077  0.3194  0.0109  0.0114  0.0095  0.0724  0.0133  0.0199  0.0312  0.4411  0.0767  0.3687  0.2369  0.1455  0.1530  0.1660  0.2169  0.3878  0.2466  0.4716  
iter 4: 0.0060  0.0018  0.7874  0.1604  0.0832  0.7394  0.0591  0.0940  0.0479  0.4926  0.0301  0.3317  0.6896  0.1434  0.0863  0.0601  0.2562  0.1227  0.0711  0.1900  0.3089  0.3194  0.0105  0.0112  0.0095  0.0707  0.0140  0.0202  0.0318  0.4447  0.0784  0.3674  0.2339  0.1423  0.1592  0.1700  0.2254  0.3886  0.2552  0.4662  

seasonalData_lag_wr_all_imp

missRanger object. Extract imputed data via $data
- best iteration: 3 
- best average OOB imputation error: 0.203846 

data_all_wr <- seasonalData_lag_wr_all_imp$data
data_all_wr$fantasyPointsMC_lag <- scale(data_all_wr$fantasyPoints_lag, scale = FALSE) # mean-centered
data_all_wr_matrix <- data_all_wr %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()
newData_wr <- data_all_wr %>% 
  filter(season == max(season, na.rm = TRUE)) %>% 
  select(-fantasyPoints_lag, -fantasyPointsMC_lag)
newData_wr_matrix <- data_all_wr_matrix[
  data_all_wr_matrix[, "season"] == max(data_all_wr_matrix[, "season"], na.rm = TRUE), # keep only rows with the most recent season
  , # all columns
  drop = FALSE]

dropCol_wr <- which(colnames(newData_wr_matrix) %in% c("fantasyPoints_lag","fantasyPointsMC_lag"))
newData_wr_matrix <- newData_wr_matrix[, -dropCol_wr, drop = FALSE]

seasonalData_lag_wr_train_imp <- missRanger::missRanger(
  seasonalData_lag_wr_train,
  pmm.k = 5,
  verbose = 2,
  seed = 52242,
  keep_forests = TRUE)

Variables to impute:        fantasy_points, fantasy_points_ppr, special_teams_tds, years_of_experience, receiving_epa, racr, air_yards_share, target_share, wopr, fantasyPoints_lag, rookie_year, rushing_epa, draft_number, gs, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, adot.rec, rat.rec, drop_percent.rec, rec_br.rec, ybc_r.rec, yac_r.rec, ybc_att.rush, yac_att.rush
Variables used to impute:   gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic, height, weight, rookie_year, draft_number, fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_epa, rushing_2pt_conversions, receptions, targets, receiving_yards, receiving_tds, receiving_fumbles, receiving_fumbles_lost, receiving_air_yards, receiving_yards_after_catch, receiving_first_downs, receiving_epa, receiving_2pt_conversions, racr, target_share, air_yards_share, wopr, special_teams_tds, ybc_att.rush, yac_att.rush, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, ybc_r.rec, yac_r.rec, adot.rec, rat.rec, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, drop_percent.rec, rec_br.rec

    fntsy_  fnts__  spcl__  yrs_f_  rcvng_  racr    ar_yr_  trgt_s  wopr    fntsP_  rok_yr  rshng_  drft_n  gs  att.rs  yds.rs  td.rsh  x1d.rs  ybc.rs  yc.rsh  brk_tkl.rs  att_b.  tgt.rc  rec.rc  yds.rc  td.rec  x1d.rc  ybc.rc  yac.rc  brk_tkl.rc  drp.rc  int.rc  adt.rc  rat.rc  drp_p.  rc_br.  ybc_r.  yc_r.r  ybc_t.  yc_tt.
iter 1: 0.0064  0.0010  0.7029  0.1611  0.1089  0.8443  0.1021  0.1643  0.0427  0.4935  0.0173  0.3461  0.6788  0.1427  0.1364  0.0993  0.2403  0.1243  0.0979  0.2745  0.6190  0.3171  0.0201  0.0147  0.0159  0.0734  0.0140  0.0280  0.0454  0.4502  0.4439  0.4733  0.3088  0.1641  0.4547  0.2192  0.2439  0.4227  0.2921  0.5068  
iter 2: 0.0063  0.0020  0.7835  0.1630  0.0901  0.8044  0.0674  0.0936  0.0479  0.4930  0.0331  0.3235  0.7090  0.1417  0.0896  0.0659  0.2659  0.1273  0.0752  0.1920  0.3068  0.3225  0.0112  0.0116  0.0101  0.0753  0.0141  0.0210  0.0333  0.4431  0.0809  0.3676  0.2571  0.1617  0.1735  0.1797  0.2441  0.3996  0.2911  0.4923  
iter 3: 0.0063  0.0020  0.7710  0.1639  0.0881  0.7954  0.0646  0.0982  0.0515  0.4956  0.0338  0.3200  0.7088  0.1413  0.0900  0.0640  0.2565  0.1250  0.0735  0.1989  0.3082  0.3280  0.0114  0.0119  0.0096  0.0763  0.0141  0.0216  0.0326  0.4388  0.0807  0.3703  0.2582  0.1623  0.1657  0.2018  0.2375  0.4016  0.2838  0.4794  
iter 4: 0.0062  0.0020  0.7792  0.1625  0.0877  0.8043  0.0632  0.0919  0.0477  0.4963  0.0341  0.3239  0.7038  0.1420  0.0950  0.0653  0.2664  0.1309  0.0767  0.2025  0.2933  0.3076  0.0109  0.0119  0.0097  0.0745  0.0143  0.0217  0.0326  0.4432  0.0804  0.3688  0.2584  0.1605  0.1931  0.2013  0.2378  0.4119  0.2824  0.4860  

seasonalData_lag_wr_train_imp

missRanger object. Extract imputed data via $data
- best iteration: 3 
- best average OOB imputation error: 0.2110456 

data_train_wr <- seasonalData_lag_wr_train_imp$data
data_train_wr$fantasyPointsMC_lag <- scale(data_train_wr$fantasyPoints_lag, scale = FALSE) # mean-centered
data_train_wr_matrix <- data_train_wr %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()

seasonalData_lag_wr_test_imp <- predict(
  object = seasonalData_lag_wr_train_imp,
  newdata = seasonalData_lag_wr_test,
  seed = 52242)

data_test_wr <- seasonalData_lag_wr_test_imp
data_test_wr_matrix <- data_test_wr %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()

# TEs
seasonalData_lag_te_all_imp <- missRanger::missRanger(
  seasonalData_lag_te_all,
  pmm.k = 5,
  verbose = 2,
  seed = 52242,
  keep_forests = TRUE)

Variables to impute:        games, ageCentered20, ageCentered20Quadratic, fantasy_points, fantasy_points_ppr, fantasyPoints, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_2pt_conversions, receptions, targets, receiving_yards, receiving_tds, receiving_fumbles, receiving_fumbles_lost, receiving_air_yards, receiving_yards_after_catch, receiving_first_downs, receiving_2pt_conversions, special_teams_tds, years_of_experience, receiving_epa, racr, air_yards_share, target_share, wopr, fantasyPoints_lag, rookie_year, draft_number, gs, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, adot.rec, rat.rec, drop_percent.rec, rec_br.rec, ybc_r.rec, yac_r.rec, rushing_epa, ybc_att.rush, yac_att.rush
Variables used to impute:   gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic, height, weight, rookie_year, draft_number, fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_epa, rushing_2pt_conversions, receptions, targets, receiving_yards, receiving_tds, receiving_fumbles, receiving_fumbles_lost, receiving_air_yards, receiving_yards_after_catch, receiving_first_downs, receiving_epa, receiving_2pt_conversions, racr, target_share, air_yards_share, wopr, special_teams_tds, ybc_att.rush, yac_att.rush, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, ybc_r.rec, yac_r.rec, adot.rec, rat.rec, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, drop_percent.rec, rec_br.rec

    games   agCn20  agC20Q  fntsy_  fnts__  fntsyP  carris  rshng_y rshng_t rshng_f rshng_fm_   rshng_fr_   rsh_2_  rcptns  targts  rcvng_y rcvng_t rcvng_f rcvng_fm_   rcvng_r_    rcv___  rcvng_fr_   rcv_2_  spcl__  yrs_f_  rcvng_p racr    ar_yr_  trgt_s  wopr    fntsP_  rok_yr  drft_n  gs  att.rs  yds.rs  td.rsh  x1d.rs  ybc.rs  yc.rsh  brk_tkl.rs  att_b.  tgt.rc  rec.rc  yds.rc  td.rec  x1d.rc  ybc.rc  yac.rc  brk_tkl.rc  drp.rc  int.rc  adt.rc  rat.rc  drp_p.  rc_br.  ybc_r.  yc_r.r  rshng_p ybc_t.  yc_tt.
iter 1: 0.8157  0.0061  0.0030  0.3406  0.0194  0.0039  0.5253  0.2259  0.2452  0.7083  0.6874  0.0802  1.1303  0.0281  0.0558  0.0255  0.0845  0.8134  0.4317  0.0655  0.0784  0.0253  0.9716  1.0271  0.1530  0.1689  0.6899  0.1092  0.4432  0.1004  0.4764  0.0180  0.6054  0.3846  0.0762  0.0832  0.1564  0.0618  0.0704  0.2123  0.3921  0.6733  0.0290  0.0207  0.0226  0.1012  0.0212  0.0420  0.0603  0.4332  0.4640  0.4996  0.2804  0.1667  0.3542  0.1652  0.2843  0.3948  0.3270  0.6620  0.7439  
iter 2: 0.1712  0.0175  0.0256  0.0106  0.0037  0.0055  0.1140  0.1113  0.0990  0.5369  0.7422  0.0852  1.1286  0.0193  0.0200  0.0128  0.0862  0.4248  0.4659  0.0206  0.0529  0.0217  0.9711  1.0114  0.1561  0.1397  0.6715  0.0766  0.1819  0.1085  0.4649  0.0366  0.6346  0.3880  0.0722  0.0728  0.1592  0.0680  0.0759  0.2034  0.3651  0.6811  0.0164  0.0158  0.0161  0.1080  0.0211  0.0327  0.0475  0.4342  0.1149  0.4173  0.2589  0.1742  0.1467  0.1531  0.2941  0.3851  0.3357  0.6846  0.7397  
iter 3: 0.1689  0.0170  0.0261  0.0114  0.0040  0.0056  0.1190  0.1155  0.0978  0.6088  0.7899  0.0945  1.1731  0.0195  0.0203  0.0132  0.0964  0.4270  0.4608  0.0202  0.0525  0.0214  0.9694  1.0265  0.1560  0.1380  0.6453  0.0751  0.1794  0.1204  0.4642  0.0364  0.6369  0.3853  0.0779  0.0786  0.1497  0.0569  0.0932  0.2027  0.4003  0.6633  0.0171  0.0164  0.0167  0.1049  0.0220  0.0335  0.0466  0.4371  0.1141  0.4304  0.2665  0.1775  0.1464  0.1537  0.2916  0.3720  0.3137  0.6640  0.7770  

seasonalData_lag_te_all_imp

missRanger object. Extract imputed data via $data
- best iteration: 2 
- best average OOB imputation error: 0.2477098 

data_all_te <- seasonalData_lag_te_all_imp$data
data_all_te$fantasyPointsMC_lag <- scale(data_all_te$fantasyPoints_lag, scale = FALSE) # mean-centered
data_all_te_matrix <- data_all_te %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()
newData_te <- data_all_te %>% 
  filter(season == max(season, na.rm = TRUE)) %>% 
  select(-fantasyPoints_lag, -fantasyPointsMC_lag)
newData_te_matrix <- data_all_te_matrix[
  data_all_te_matrix[, "season"] == max(data_all_te_matrix[, "season"], na.rm = TRUE), # keep only rows with the most recent season
  , # all columns
  drop = FALSE]

dropCol_te <- which(colnames(newData_te_matrix) %in% c("fantasyPoints_lag","fantasyPointsMC_lag"))
newData_te_matrix <- newData_te_matrix[, -dropCol_te, drop = FALSE]

seasonalData_lag_te_train_imp <- missRanger::missRanger(
  seasonalData_lag_te_train,
  pmm.k = 5,
  verbose = 2,
  seed = 52242,
  keep_forests = TRUE)

Variables to impute:        games, years_of_experience, ageCentered20, ageCentered20Quadratic, fantasy_points, fantasy_points_ppr, fantasyPoints, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_2pt_conversions, receptions, targets, receiving_yards, receiving_tds, receiving_fumbles, receiving_fumbles_lost, receiving_air_yards, receiving_yards_after_catch, receiving_first_downs, receiving_2pt_conversions, special_teams_tds, receiving_epa, racr, air_yards_share, target_share, wopr, fantasyPoints_lag, rookie_year, draft_number, gs, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, adot.rec, rat.rec, drop_percent.rec, rec_br.rec, ybc_r.rec, yac_r.rec, rushing_epa, ybc_att.rush, yac_att.rush
Variables used to impute:   gsis_id, season, games, gs, years_of_experience, age, ageCentered20, ageCentered20Quadratic, height, weight, rookie_year, draft_number, fantasy_points, fantasy_points_ppr, fantasyPoints, fantasyPoints_lag, carries, rushing_yards, rushing_tds, rushing_fumbles, rushing_fumbles_lost, rushing_first_downs, rushing_epa, rushing_2pt_conversions, receptions, targets, receiving_yards, receiving_tds, receiving_fumbles, receiving_fumbles_lost, receiving_air_yards, receiving_yards_after_catch, receiving_first_downs, receiving_epa, receiving_2pt_conversions, racr, target_share, air_yards_share, wopr, special_teams_tds, ybc_att.rush, yac_att.rush, att.rush, yds.rush, td.rush, x1d.rush, ybc.rush, yac.rush, brk_tkl.rush, att_br.rush, ybc_r.rec, yac_r.rec, adot.rec, rat.rec, tgt.rec, rec.rec, yds.rec, td.rec, x1d.rec, ybc.rec, yac.rec, brk_tkl.rec, drop.rec, int.rec, drop_percent.rec, rec_br.rec

    games   yrs_f_  agCn20  agC20Q  fntsy_  fnts__  fntsyP  carris  rshng_y rshng_t rshng_f rshng_fm_   rshng_fr_   rsh_2_  rcptns  targts  rcvng_y rcvng_t rcvng_f rcvng_fm_   rcvng_r_    rcv___  rcvng_fr_   rcv_2_  spcl__  rcvng_p racr    ar_yr_  trgt_s  wopr    fntsP_  rok_yr  drft_n  gs  att.rs  yds.rs  td.rsh  x1d.rs  ybc.rs  yc.rsh  brk_tkl.rs  att_b.  tgt.rc  rec.rc  yds.rc  td.rec  x1d.rc  ybc.rc  yac.rc  brk_tkl.rc  drp.rc  int.rc  adt.rc  rat.rc  drp_p.  rc_br.  ybc_r.  yc_r.r  rshng_p ybc_t.  yc_tt.
iter 1: 0.8094  0.1093  0.0070  0.0035  0.3272  0.0235  0.0052  0.2840  0.1426  0.2634  0.8628  0.7885  0.0924  1.1067  0.0298  0.0611  0.0249  0.0969  0.8177  0.4537  0.0650  0.0804  0.0249  0.9680  1.0235  0.1738  0.5438  0.0868  0.4123  0.1172  0.4597  0.0189  0.6057  0.3973  0.0877  0.0886  0.1516  0.0467  0.0593  0.2086  0.4018  0.6464  0.0296  0.0223  0.0237  0.1062  0.0207  0.0428  0.0579  0.4367  0.4700  0.4818  0.3045  0.1724  0.4722  0.2410  0.2693  0.4025  0.3943  0.4791  0.7521  
iter 2: 0.1728  0.1469  0.0179  0.0289  0.0104  0.0039  0.0051  0.0863  0.0763  0.1528  0.7880  0.9474  0.0849  1.0234  0.0193  0.0198  0.0137  0.0915  0.4327  0.4835  0.0238  0.0558  0.0221  0.9617  1.0376  0.1464  0.5141  0.0630  0.1827  0.1062  0.4562  0.0379  0.6361  0.3908  0.0641  0.0767  0.1425  0.0603  0.0747  0.1970  0.3903  0.6647  0.0182  0.0171  0.0165  0.1074  0.0226  0.0332  0.0503  0.4361  0.1170  0.4096  0.2673  0.1862  0.2255  0.2386  0.2793  0.4004  0.3648  0.5070  0.7621  
iter 3: 0.1713  0.1447  0.0195  0.0276  0.0104  0.0036  0.0051  0.0796  0.0889  0.1611  0.8505  0.9348  0.0904  1.0447  0.0196  0.0205  0.0134  0.0901  0.4465  0.4867  0.0233  0.0569  0.0222  0.9519  1.0103  0.1457  0.5062  0.0617  0.1698  0.1115  0.4530  0.0382  0.6521  0.3899  0.0665  0.0681  0.1457  0.0647  0.0866  0.2055  0.3919  0.6791  0.0169  0.0169  0.0168  0.1107  0.0213  0.0339  0.0500  0.4315  0.1200  0.4148  0.2745  0.1822  0.1947  0.2205  0.2778  0.4032  0.3639  0.4933  0.7741  

seasonalData_lag_te_train_imp

missRanger object. Extract imputed data via $data
- best iteration: 2 
- best average OOB imputation error: 0.2519559 

data_train_te <- seasonalData_lag_te_train_imp$data
data_train_te$fantasyPointsMC_lag <- scale(data_train_te$fantasyPoints_lag, scale = FALSE) # mean-centered
data_train_te_matrix <- data_train_te %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()

seasonalData_lag_te_test_imp <- predict(
  object = seasonalData_lag_te_train_imp,
  newdata = seasonalData_lag_te_test,
  seed = 52242)

data_test_te <- seasonalData_lag_te_test_imp
data_test_te_matrix <- data_test_te %>%
  mutate(across(where(is.factor), ~ as.numeric(as.integer(.)))) %>% 
  as.matrix()

19.5 Identify Cores for Parallel Processing

num_cores <- parallel::detectCores() - 1
num_true_cores <- parallel::detectCores(logical = FALSE) - 1

num_cores

[1] 4

We use the future (Bengtsson, 2025) package for parallel (faster) processing.

future::plan(future::multisession, workers = num_cores)

19.6 Fitting the Traditional Regression Models

19.6.1 Regression with One Predictor

# Set seed for reproducibility
set.seed(52242)

# Set up Cross-Validation
folds <- rsample::group_vfold_cv(
  data_train_qb,
  group = gsis_id, # ensures all rows for a player are in the training set or all in the validation set for each fold
  v = 10) # 10-fold cross-validation

# Define Recipe (Formula)
rec <- recipes::recipe(
  fantasyPoints_lag ~ ageCentered20,
  data = data_train_qb)

# Define Model
lm_spec <- parsnip::linear_reg() %>%
  parsnip::set_engine("lm") %>%
  parsnip::set_mode("regression")

# Workflow
lm_wf <- workflows::workflow() %>%
  workflows::add_recipe(rec) %>%
  workflows::add_model(lm_spec)

# Fit Model with Cross-Validation
cv_results <- tune::fit_resamples(
  lm_wf,
  resamples = folds,
  metrics = yardstick::metric_set(rmse, mae, rsq),
  control = tune::control_resamples(save_pred = TRUE)
)

# View Cross-Validation metrics
tune::collect_metrics(cv_results)

# Fit Final Model on Training Data
final_model <- workflows::fit(
  lm_wf,
  data = data_train_qb)

# View Coefficients
final_model %>% 
  workflows::extract_fit_parsnip() %>% 
  broom::tidy()

# Predict on Test Data
df <- data_test_qb %>%
  mutate(pred = predict(final_model, new_data = data_test_qb)$.pred)

# Evaluate Accuracy of Predictions
petersenlab::accuracyOverall(
  predicted = df$pred,
  actual = df$fantasyPoints_lag,
  dropUndefined = TRUE
)

# Calculate combined range for axes
axis_limits <- range(c(df$pred, df$fantasyPoints_lag), na.rm = TRUE)

ggplot(
  df,
  aes(
    x = pred,
    y = fantasyPoints_lag)) +
  geom_point(
    size = 2,
    alpha = 0.6) +
  geom_abline(
    slope = 1,
    intercept = 0,
    color = "blue",
    linetype = "dashed") +
  coord_equal(
    xlim = axis_limits,
    ylim = axis_limits) +
  labs(
    title = "Predicted vs Actual Fantasy Points (Test Data)",
    x = "Predicted Fantasy Points",
    y = "Actual Fantasy Points"
  ) +
  theme_classic()

Figure 19.1: Predicted Versus Actual Fantasy Points for Regression Model with One Predictor (Player Age).

newData_qb %>%
  mutate(fantasyPoints_lag = predict(final_model, new_data = newData_qb)$.pred) %>% 
  left_join(
    .,
    nfl_playerIDs %>% select(gsis_id, name),
    by = "gsis_id"
  ) %>% 
  select(name, fantasyPoints_lag) %>% 
  arrange(-fantasyPoints_lag)

19.6.2 Regression with Multiple Predictors

# Set seed for reproducibility
set.seed(52242)

# Set up Cross-Validation
folds <- rsample::group_vfold_cv(
  data_train_qb,
  group = gsis_id, # ensures all rows for a player are in the training set or all in the validation set for each fold
  v = 10) # 10-fold cross-validation

# Define Recipe (Formula)
rec <- recipes::recipe(
  fantasyPoints_lag ~ ., # use all predictors
  data = data_train_qb %>% select(-gsis_id, -fantasyPointsMC_lag))

# Define Model
lm_spec <- parsnip::linear_reg() %>%
  parsnip::set_engine("lm") %>%
  parsnip::set_mode("regression")

# Workflow
lm_wf <- workflows::workflow() %>%
  workflows::add_recipe(rec) %>%
  workflows::add_model(lm_spec)

# Fit Model with Cross-Validation
cv_results <- tune::fit_resamples(
  lm_wf,
  resamples = folds,
  metrics = yardstick::metric_set(rmse, mae, rsq),
  control = tune::control_resamples(save_pred = TRUE)
)

# View Cross-Validation metrics
tune::collect_metrics(cv_results)

# Fit Final Model on Training Data
final_model <- workflows::fit(
  lm_wf,
  data = data_train_qb)

# View Coefficients
final_model %>% 
  workflows::extract_fit_parsnip() %>% 
  broom::tidy()

# Predict on Test Data
df <- data_test_qb %>%
  mutate(pred = predict(final_model, new_data = data_test_qb)$.pred)

# Evaluate Accuracy of Predictions
petersenlab::accuracyOverall(
  predicted = df$pred,
  actual = df$fantasyPoints_lag,
  dropUndefined = TRUE
)

# Calculate combined range for axes
axis_limits <- range(c(df$pred, df$fantasyPoints_lag), na.rm = TRUE)

ggplot(
  df,
  aes(
    x = pred,
    y = fantasyPoints_lag)) +
  geom_point(
    size = 2,
    alpha = 0.6) +
  geom_abline(
    slope = 1,
    intercept = 0,
    color = "blue",
    linetype = "dashed") +
  coord_equal(
    xlim = axis_limits,
    ylim = axis_limits) +
  labs(
    title = "Predicted vs Actual Fantasy Points (Test Data)",
    x = "Predicted Fantasy Points",
    y = "Actual Fantasy Points"
  ) +
  theme_classic()

Figure 19.2: Predicted Versus Actual Fantasy Points for Regression Model with Multiple Predictors.

newData_qb %>%
  mutate(fantasyPoints_lag = predict(final_model, new_data = newData_qb)$.pred) %>% 
  left_join(
    .,
    nfl_playerIDs %>% select(gsis_id, name),
    by = "gsis_id"
  ) %>% 
  select(name, fantasyPoints_lag) %>% 
  arrange(-fantasyPoints_lag)

19.7 Fitting the Machine Learning Models

19.7.1 Least Absolute Shrinkage and Selection Option (LASSO)

# Set seed for reproducibility
set.seed(52242)

# Set up Cross-Validation
folds <- rsample::group_vfold_cv(
  data_train_qb,
  group = gsis_id, # ensures all rows for a player are in the training set or all in the validation set for each fold
  v = 10) # 10-fold cross-validation

# Define Recipe (Formula)
rec <- recipes::recipe(
  fantasyPoints_lag ~ ., # use all predictors
  data = data_train_qb %>% select(-gsis_id, -fantasyPointsMC_lag))

# Define Model
lasso_spec <- 
  parsnip::linear_reg(
    penalty = tune::tune(),
    mixture = 1) %>%
  parsnip::set_engine("glmnet")

# Workflow
lasso_wf <- workflows::workflow() %>%
  workflows::add_recipe(rec) %>%
  workflows::add_model(lasso_spec)

# Define grid of penalties to try (log scale is typical)
penalty_grid <- dials::grid_regular(
  dials::penalty(range = c(-4, -1)),
  levels = 20)

# Tune the Penalty Parameter
cv_results <- tune::tune_grid(
  lasso_wf,
  resamples = folds,
  grid = penalty_grid,
  metrics = yardstick::metric_set(rmse, mae, rsq),
  control = tune::control_grid(save_pred = TRUE)
)

# View Cross-Validation metrics
tune::collect_metrics(cv_results)

# Identify best penalty
tune::select_best(cv_results, metric = "rmse")

tune::select_best(cv_results, metric = "mae")

tune::select_best(cv_results, metric = "rsq")

best_penalty <- tune::select_best(cv_results, metric = "mae")

# Finalize Workflow with Best Penalty
final_wf <- tune::finalize_workflow(
  lasso_wf,
  best_penalty)

# Fit Final Model on Training Data
final_model <- workflows::fit(
  final_wf,
  data = data_train_qb)

# View Coefficients
final_model %>% 
  workflows::extract_fit_parsnip() %>% 
  broom::tidy()

# Predict on Test Data
df <- data_test_qb %>%
  mutate(pred = predict(final_model, new_data = data_test_qb)$.pred)

# Evaluate Accuracy of Predictions
petersenlab::accuracyOverall(
  predicted = df$pred,
  actual = df$fantasyPoints_lag,
  dropUndefined = TRUE
)

# Calculate combined range for axes
axis_limits <- range(c(df$pred, df$fantasyPoints_lag), na.rm = TRUE)

ggplot(
  df,
  aes(
    x = pred,
    y = fantasyPoints_lag)) +
  geom_point(
    size = 2,
    alpha = 0.6) +
  geom_abline(
    slope = 1,
    intercept = 0,
    color = "blue",
    linetype = "dashed") +
  coord_equal(
    xlim = axis_limits,
    ylim = axis_limits) +
  labs(
    title = "Predicted vs Actual Fantasy Points (Test Data)",
    x = "Predicted Fantasy Points",
    y = "Actual Fantasy Points"
  ) +
  theme_classic()

Figure 19.3: Predicted Versus Actual Fantasy Points for Least Absolute Shrinkage and Selection Option (LASSO) Model.

newData_qb %>%
  mutate(fantasyPoints_lag = predict(final_model, new_data = newData_qb)$.pred) %>% 
  left_join(
    .,
    nfl_playerIDs %>% select(gsis_id, name),
    by = "gsis_id"
  ) %>% 
  select(name, fantasyPoints_lag) %>% 
  arrange(-fantasyPoints_lag)

19.7.2 Ridge Regression

# Set seed for reproducibility
set.seed(52242)

# Set up Cross-Validation
folds <- rsample::group_vfold_cv(
  data_train_qb,
  group = gsis_id, # ensures all rows for a player are in the training set or all in the validation set for each fold
  v = 10) # 10-fold cross-validation

# Define Recipe (Formula)
rec <- recipes::recipe(
  fantasyPoints_lag ~ ., # use all predictors
  data = data_train_qb %>% select(-gsis_id, -fantasyPointsMC_lag))

# Define Model
ridge_spec <- 
  parnsnip::linear_reg(
    penalty = tune::tune(),
    mixture = 0) %>%
  parsnip::set_engine("glmnet")

Error in loadNamespace(x): there is no package called 'parnsnip'

# Workflow
ridge_wf <- workflows::workflow() %>%
  workflows::add_recipe(rec) %>%
  workflows::add_model(ridge_spec)

Error: object 'ridge_spec' not found

# Define grid of penalties to try (log scale is typical)
penalty_grid <- dials::grid_regular(
  dials::penalty(range = c(-4, -1)),
  levels = 20)

# Tune the Penalty Parameter
cv_results <- tune::tune_grid(
  ridge_wf,
  resamples = folds,
  grid = penalty_grid,
  metrics = yardstick::metric_set(rmse, mae, rsq),
  control = tune::control_grid(save_pred = TRUE)
)

Error: object 'ridge_wf' not found

# View Cross-Validation metrics
tune::collect_metrics(cv_results)

# Identify best penalty
tune::select_best(cv_results, metric = "rmse")

tune::select_best(cv_results, metric = "mae")

tune::select_best(cv_results, metric = "rsq")

best_penalty <- tune::select_best(cv_results, metric = "mae")

# Finalize Workflow with Best Penalty
final_wf <- tune::finalize_workflow(
  ridge_wf,
  best_penalty)

Error: object 'ridge_wf' not found

# Fit Final Model on Training Data
final_model <- workflows::fit(
  final_wf,
  data = data_train_qb)

# View Coefficients
final_model %>% 
  workflows::extract_fit_parsnip() %>% 
  broom::tidy()

# Predict on Test Data
df <- data_test_qb %>%
  mutate(pred = predict(final_model, new_data = data_test_qb)$.pred)

# Evaluate Accuracy of Predictions
petersenlab::accuracyOverall(
  predicted = df$pred,
  actual = df$fantasyPoints_lag,
  dropUndefined = TRUE
)

# Calculate combined range for axes
axis_limits <- range(c(df$pred, df$fantasyPoints_lag), na.rm = TRUE)

ggplot(
  df,
  aes(
    x = pred,
    y = fantasyPoints_lag)) +
  geom_point(
    size = 2,
    alpha = 0.6) +
  geom_abline(
    slope = 1,
    intercept = 0,
    color = "blue",
    linetype = "dashed") +
  coord_equal(
    xlim = axis_limits,
    ylim = axis_limits) +
  labs(
    title = "Predicted vs Actual Fantasy Points (Test Data)",
    x = "Predicted Fantasy Points",
    y = "Actual Fantasy Points"
  ) +
  theme_classic()

Figure 19.4: Predicted Versus Actual Fantasy Points for Ridge Regression Model.

newData_qb %>%
  mutate(fantasyPoints_lag = predict(final_model, new_data = newData_qb)$.pred) %>% 
  left_join(
    .,
    nfl_playerIDs %>% select(gsis_id, name),
    by = "gsis_id"
  ) %>% 
  select(name, fantasyPoints_lag) %>% 
  arrange(-fantasyPoints_lag)

19.7.3 Elastic Net

# Set seed for reproducibility
set.seed(52242)

# Set up Cross-Validation
folds <- rsample::group_vfold_cv(
  data_train_qb,
  group = gsis_id, # ensures all rows for a player are in the training set or all in the validation set for each fold
  v = 10) # 10-fold cross-validation

# Define Recipe (Formula)
rec <- recipes::recipe(
  fantasyPoints_lag ~ ., # use all predictors
  data = data_train_qb %>% select(-gsis_id, -fantasyPointsMC_lag))

# Define Model
enet_spec <- 
  parsnip::linear_reg(
    penalty = tune::tune(),
    mixture = tune::tune()) %>%
  parsnip::set_engine("glmnet")

# Workflow
enet_wf <- workflows::workflow() %>%
  workflows::add_recipe(rec) %>%
  workflows::add_model(enet_spec)

# Define a regular grid for both penalty and mixture
grid_enet <- dials::grid_regular(
  dials::penalty(range = c(-4, -1)),
  dials::mixture(range = c(0, 1)),
  levels = c(20, 5) # 20 penalty values × 5 mixture values
)

# Tune the Grid
cv_results <- tune::tune_grid(
  enet_wf,
  resamples = folds,
  grid = grid_enet,
  metrics = yardstick::metric_set(rmse, mae, rsq),
  control = tune::control_grid(save_pred = TRUE)
)

# View Cross-Validation metrics
tune::collect_metrics(cv_results)

# Identify best penalty
tune::select_best(cv_results, metric = "rmse")

tune::select_best(cv_results, metric = "mae")

tune::select_best(cv_results, metric = "rsq")

best_penalty <- tune::select_best(cv_results, metric = "mae")

# Finalize Workflow with Best Penalty
final_wf <- tune::finalize_workflow(
  enet_wf,
  best_penalty)

# Fit Final Model on Training Data
final_model <- workflows::fit(
  final_wf,
  data = data_train_qb)

# View Coefficients
final_model %>% 
  workflows::extract_fit_parsnip() %>% 
  broom::tidy()

# Predict on Test Data
df <- data_test_qb %>%
  mutate(pred = predict(final_model, new_data = data_test_qb)$.pred)

# Evaluate Accuracy of Predictions
petersenlab::accuracyOverall(
  predicted = df$pred,
  actual = df$fantasyPoints_lag,
  dropUndefined = TRUE
)

# Calculate combined range for axes
axis_limits <- range(c(df$pred, df$fantasyPoints_lag), na.rm = TRUE)

ggplot(
  df,
  aes(
    x = pred,
    y = fantasyPoints_lag)) +
  geom_point(
    size = 2,
    alpha = 0.6) +
  geom_abline(
    slope = 1,
    intercept = 0,
    color = "blue",
    linetype = "dashed") +
  coord_equal(
    xlim = axis_limits,
    ylim = axis_limits) +
  labs(
    title = "Predicted vs Actual Fantasy Points (Test Data)",
    x = "Predicted Fantasy Points",
    y = "Actual Fantasy Points"
  ) +
  theme_classic()

Figure 19.5: Predicted Versus Actual Fantasy Points for Elastic Net Model.

newData_qb %>%
  mutate(fantasyPoints_lag = predict(final_model, new_data = newData_qb)$.pred) %>% 
  left_join(
    .,
    nfl_playerIDs %>% select(gsis_id, name),
    by = "gsis_id"
  ) %>% 
  select(name, fantasyPoints_lag) %>% 
  arrange(-fantasyPoints_lag)

19.7.4 Random Forest Machine Learning

19.7.4.1 Cross-Sectional Data

# Set seed for reproducibility
set.seed(52242)

# Set up Cross-Validation
folds <- rsample::group_vfold_cv(
  data_train_qb,
  group = gsis_id, # ensures all rows for a player are in the training set or all in the validation set for each fold
  v = 10) # 10-fold cross-validation

# Define Recipe (Formula)
rec <- recipes::recipe(
  fantasyPoints_lag ~ ., # use all predictors
  data = data_train_qb %>% select(-gsis_id, -fantasyPointsMC_lag))

# Define Model
rf_spec <- 
  parsnip::rand_forest(
    mtry = tune::tune(),
    min_n = tune::tune(),
    trees = 500) %>%
  parsnip::set_mode("regression") %>%
  parsnip::set_engine(
    "ranger",
    importance = "impurity")

# Workflow
rf_wf <- workflows::workflow() %>%
  workflows::add_recipe(rec) %>%
  workflows::add_model(rf_spec)

# Create Grid
n_predictors <- recipes::prep(rec) %>%
  recipes::juice() %>%
  dplyr::select(-fantasyPoints_lag) %>%
  ncol()

# Dynamically define ranges based on data
rf_params <- hardhat::extract_parameter_set_dials(rf_spec) %>%
  dials:::update.parameters(
    mtry = dials::mtry(range = c(1L, n_predictors)),
    min_n = dials::min_n(range = c(2L, 10L))
  )

rf_grid <- dials::grid_random(rf_params, size = 15) #dials::grid_regular(rf_params, levels = 5)

# Tune the Grid
cv_results <- tune::tune_grid(
  rf_wf,
  resamples = folds,
  grid = rf_grid,
  metrics = yardstick::metric_set(rmse, mae, rsq),
  control = tune::control_grid(save_pred = TRUE)
)

# View Cross-Validation metrics
tune::collect_metrics(cv_results)

# Identify best penalty
tune::select_best(cv_results, metric = "rmse")

tune::select_best(cv_results, metric = "mae")

tune::select_best(cv_results, metric = "rsq")

best_penalty <- tune::select_best(cv_results, metric = "mae")

# Finalize Workflow with Best Penalty
final_wf <- tune::finalize_workflow(
  rf_wf,
  best_penalty)

# Fit Final Model on Training Data
final_model <- workflows::fit(
  final_wf,
  data = data_train_qb)

# View Feature Importance
rf_fit <- final_model %>% 
  workflows::extract_fit_parsnip()

rf_fit

parsnip model object

Ranger result

Call:
 ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~3L,      x), num.trees = ~500, min.node.size = min_rows(~9L, x), importance = ~"impurity",      num.threads = 1, verbose = FALSE, seed = sample.int(10^5,          1)) 

Type:                             Regression 
Number of trees:                  500 
Sample size:                      1582 
Number of independent variables:  73 
Mtry:                             3 
Target node size:                 9 
Variable importance mode:         impurity 
Splitrule:                        variance 
OOB prediction error (MSE):       6315.315 
R squared (OOB):                  0.5148104 

ranger_obj <- rf_fit$fit

ranger_obj$variable.importance

                   season                     games                        gs 
              126359.5422               319578.3505               382535.6517 
      years_of_experience                       age             ageCentered20 
              124783.8817               199342.1467               207939.2802 
   ageCentered20Quadratic                    height                    weight 
              205495.8993                73752.5223               135823.1306 
              rookie_year              draft_number            fantasy_points 
              121270.5272               235944.0545               641663.0852 
       fantasy_points_ppr             fantasyPoints               completions 
              753816.3235               673276.9667               577883.0044 
                 attempts             passing_yards               passing_tds 
              469528.4582               692904.1588               535369.7048 
    passing_interceptions            sacks_suffered           sack_yards_lost 
              154679.5322               294698.8048               237255.5206 
             sack_fumbles         sack_fumbles_lost         passing_air_yards 
               92927.1076                52731.7059               389664.7190 
passing_yards_after_catch       passing_first_downs               passing_epa 
              371266.0798               533911.4945               469375.8940 
             passing_cpoe   passing_2pt_conversions                      pacr 
              159111.5776                43285.8990               136591.0157 
                  carries             rushing_yards               rushing_tds 
              191763.4732               227786.6905                83784.6521 
          rushing_fumbles      rushing_fumbles_lost       rushing_first_downs 
               62378.0789                31328.7581               218593.7055 
              rushing_epa   rushing_2pt_conversions         special_teams_tds 
              174856.7482                22130.6589                  432.8602 
         pocket_time.pass        pass_attempts.pass           throwaways.pass 
               88617.4566               481207.3858               330967.1992 
              spikes.pass                drops.pass           bad_throws.pass 
               74008.1352               364123.0798               411395.3664 
       times_blitzed.pass        times_hurried.pass            times_hit.pass 
              481435.4648               374116.1093               266521.6649 
     times_pressured.pass         batted_balls.pass        on_tgt_throws.pass 
              424836.7389               151484.9543               402761.4106 
           rpo_plays.pass            rpo_yards.pass         rpo_pass_att.pass 
              192637.5465               225659.8819               243476.9956 
      rpo_pass_yards.pass         rpo_rush_att.pass       rpo_rush_yards.pass 
              205064.4104                77907.9318               117492.6716 
         pa_pass_att.pass        pa_pass_yards.pass             drop_pct.pass 
              347266.3891               454836.7853               129912.7438 
       bad_throw_pct.pass           on_tgt_pct.pass         pressure_pct.pass 
              145508.9586               106607.2086               126197.3464 
             ybc_att.rush              yac_att.rush                  att.rush 
              139772.9950               120705.2307               367870.7638 
                 yds.rush                   td.rush                  x1d.rush 
              208163.2897               125274.5721               221642.2967 
                 ybc.rush                  yac.rush              brk_tkl.rush 
              188152.9223               173028.6991                52392.4196 
              att_br.rush 
               86667.3297 

# Predict on Test Data
df <- data_test_qb %>%
  mutate(pred = predict(final_model, new_data = data_test_qb)$.pred)

# Evaluate Accuracy of Predictions
petersenlab::accuracyOverall(
  predicted = df$pred,
  actual = df$fantasyPoints_lag,
  dropUndefined = TRUE
)

# Calculate combined range for axes
axis_limits <- range(c(df$pred, df$fantasyPoints_lag), na.rm = TRUE)

ggplot(
  df,
  aes(
    x = pred,
    y = fantasyPoints_lag)) +
  geom_point(
    size = 2,
    alpha = 0.6) +
  geom_abline(
    slope = 1,
    intercept = 0,
    color = "blue",
    linetype = "dashed") +
  coord_equal(
    xlim = axis_limits,
    ylim = axis_limits) +
  labs(
    title = "Predicted vs Actual Fantasy Points (Test Data)",
    x = "Predicted Fantasy Points",
    y = "Actual Fantasy Points"
  ) +
  theme_classic()

Figure 19.6: Predicted Versus Actual Fantasy Points for Random Forest Model.

newData_qb %>%
  mutate(fantasyPoints_lag = predict(final_model, new_data = newData_qb)$.pred) %>% 
  left_join(
    .,
    nfl_playerIDs %>% select(gsis_id, name),
    by = "gsis_id"
  ) %>% 
  select(name, fantasyPoints_lag) %>% 
  arrange(-fantasyPoints_lag)

Now we can stop the parallel backend:

future::plan(future::sequential)

19.7.4.2 Longitudinal Data

Approaches to estimating random forest models with longitudinal data are described by @Hu & Szymczak (2023). Below, we fit longitudinal random forest models using the MERF() function of the LongituRF package (Capitaine, 2020).

smerf <- LongituRF::MERF(
  X = data_train_qb_matrix %>% as_tibble() %>% dplyr::select(season:att_br.rush) %>% as.matrix(),
  Y = data_train_qb_matrix[,c("fantasyPoints_lag")] %>% as.matrix(),
  Z = data_train_qb_matrix[,c("pacr")] %>% as.matrix(),
  id = data_train_qb_matrix[,c("gsis_id")] %>% as.matrix(),
  time = data_train_qb_matrix[,c("ageCentered20")] %>% as.matrix(),
  ntree = 500,
  sto = "BM")

[1] "stopped after 10 iterations."

smerf$forest # the fitted random forest (obtained at the last iteration)

Call:
 randomForest(x = X, y = ystar, ntree = ntree, mtry = mtry, importance = TRUE) 
               Type of random forest: regression
                     Number of trees: 500
No. of variables tried at each split: 25

          Mean of squared residuals: 165.2308
                    % Var explained: 98.57

smerf$random_effects %>% data.frame() # the predicted random effects for each player

smerf$omega %>% data.frame() # the predicted stochastic processes

smerf$OOB # OOB error at each iteration

 [1] 178.3309 127.1757 137.5047 151.4621 150.9636 168.9314 144.5050 155.6409
 [9] 162.2556 165.2308

# Predict on Test Data
predict(
  smerf,
  X = data_test_qb %>% as_tibble() %>% dplyr::select(season:att_br.rush) %>% as.matrix(),
  Z = data_test_qb[,c("pacr")] %>% as.matrix(),
  id = data_test_qb[,c("gsis_id")] %>% as.matrix(),
  time = data_test_qb[,c("ageCentered20")] %>% as.matrix())

Error in Ypred[w] <- f[w] + Z[w, , drop = FALSE] %*% object$random_effects[k, : replacement has length zero

plot(smerf$Vraisemblance)

Figure 19.7: Evolution of the Log-Likelihood.

19.7.5 k-Fold Cross-Validation

19.7.6 Leave-One-Out (LOO) Cross-Validation

19.7.7 Combining Tree-Boosting with Mixed Models

To combine tree-boosting with mixed models, we use the gpboost package (Sigrist et al., 2025).

Adapted from here: https://towardsdatascience.com/mixed-effects-machine-learning-for-longitudinal-panel-data-with-gpboost-part-iii-523bb38effc

19.7.7.1 Process Data

If using a gamma distribution, it requires positive-only values:

data_train_qb_matrix[,"fantasyPoints_lag"][data_train_qb_matrix[,"fantasyPoints_lag"] <= 0] <- 0.01

19.7.7.2 Specify Predictor Variables

pred_vars_qb <- data_train_qb_matrix %>% 
  as_tibble() %>% 
  select(-fantasyPoints_lag, -fantasyPointsMC_lag, -ageCentered20, ageCentered20Quadratic) %>% # -gsis_id
  names()

pred_vars_qb_categorical <- "gsis_id" # to specify categorical predictors

19.7.7.3 Specify General Model Options

model_likelihood <- "gamma" # gaussian
nrounds <- 2000 # maximum number of boosting iterations (i.e., number of trees built sequentially); more rounds = potentially better learning, but also greater risk of overfitting

19.7.7.4 Identify Optimal Tuning Parameters

For identifying the optimal tuning parameters for boosting, we partition the training data into inner training data and validation data. We randomly split the training data into 80% inner training data and 20% held-out validation data. We then use the mean absolute error as our index of prediction accuracy on the held-out validation data.

# Partition training data into inner training data and validation data
ntrain_qb <- dim(data_train_qb_matrix)[1]

set.seed(52242)
valid_tune_idx_qb <- sample.int(ntrain_qb, as.integer(0.2*ntrain_qb)) # 

folds_qb <- list(valid_tune_idx_qb)

# Specify parameter grid, gp_model, and gpb.Dataset
param_grid_qb <- list(
  "learning_rate" = c(0.2, 0.1, 0.05, 0.01), # the step size used when updating predictions after each boosting round (high values make big updates, which can speed up learning but risk overshooting; low values are usually more accurate but require more rounds)
  "max_depth" = c(3, 5, 7), # maximum depth (levels) of each decision tree; deeper trees capture more complex patterns and interactions but risk overfitting; shallower trees tend to generalize better
  "min_data_in_leaf" = c(10, 50, 100), # minimum number of training examples in a leaf node; higher values = more regularization (simpler trees)
  "lambda_l2" = c(0, 1, 5)) # L2 regularization penalty for large weights in tree splits; adds a "cost" for complexity; helps prevent overfitting by shrinking the contribution of each tree

other_params_qb <- list(
  num_leaves = 2^6) # maximum number of leaves per tree; controls the maximum complexity of each tree (along with max_depth); more leaves = more expressive models, but can overfit if min_data_in_leaf is too small; num_leaves must be consistent with max_depth, because deeper trees naturally support more leaves; max is: 2^n, where n is the largest max_depth

gp_model_qb <- gpboost::GPModel(
  group_data = data_train_qb_matrix[,"gsis_id"],
  likelihood = model_likelihood,
  group_rand_coef_data = cbind(
    data_train_qb_matrix[,"ageCentered20"],
    data_train_qb_matrix[,"ageCentered20Quadratic"]),
  ind_effect_group_rand_coef = c(1,1))

gp_data_qb <- gpboost::gpb.Dataset(
  data = data_train_qb_matrix[,pred_vars_qb],
  categorical_feature = pred_vars_qb_categorical,
  label = data_train_qb_matrix[,"fantasyPoints_lag"]) # could instead use mean-centered variable (fantasyPointsMC_lag) and add mean back afterward

# Find optimal tuning parameters
opt_params_qb <- gpboost::gpb.grid.search.tune.parameters(
  param_grid = param_grid_qb,
  params = other_params_qb,
  num_try_random = NULL,
  folds = folds_qb,
  data = gp_data_qb,
  gp_model = gp_model_qb,
  nrounds = nrounds,
  early_stopping_rounds = 50, # stops training early if the model hasn’t improved on the validation set in 50 rounds; prevents overfitting and saves time
  verbose_eval = 1,
  metric = "mae")

Error in fd$booster$update(fobj = fobj): [GPBoost] [Fatal] Inf occured in gradient wrt covariance / auxiliary parameter number 3 (counting starts at 1, total nb. par. = 4) 

opt_params_qb

Error: object 'opt_params_qb' not found

A learning rate of 1 is very high for boosting. Even if a learning rate of 1 did well in tuning, I use a lower learning rate (0.1) to avoid overfitting. I also added some light regularization (lambda_l2) for better generalization. I also set the maximum tree depth (max_depth) at 5 to capture complex (up to 5-way) interactions, and set the maximum number of terminal nodes (num_leaves) per tree at 2^5 (32). I set the minimum number of samples in any leaf (min_data_in_leaf) to be 10.

19.7.7.5 Specify Model and Tuning Parameters

gp_model_qb <- gpboost::GPModel(
  group_data = data_train_qb_matrix[,"gsis_id"],
  likelihood = model_likelihood,
  group_rand_coef_data = cbind(
    data_train_qb_matrix[,"ageCentered20"],
    data_train_qb_matrix[,"ageCentered20Quadratic"]),
  ind_effect_group_rand_coef = c(1,1))

gp_data_qb <- gpboost::gpb.Dataset(
  data = data_train_qb_matrix[,pred_vars_qb],
  categorical_feature = pred_vars_qb_categorical,
  label = data_train_qb_matrix[,"fantasyPoints_lag"])

params_qb <- list(
  learning_rate = 0.1,
  max_depth = 5,
  min_data_in_leaf = 10,
  lambda_l2 = 1,
  num_leaves = 2^5,
  num_threads = num_cores)

nrounds_qb <- 123 # identify optimal number of trees through iteration and cross-validation

#gp_model_qb$set_optim_params(params = list(optimizer_cov = "nelder_mead")) # to speed up model estimation

19.7.7.6 Fit Model

gp_model_fit_qb <- gpboost::gpb.train(
  data = gp_data_qb,
  gp_model = gp_model_qb,
  nrounds = nrounds_qb,
  params = params_qb) # verbose = 0

[GPBoost] [Info] Total Bins 8709
[GPBoost] [Info] Number of data points in the train set: 1582, number of used features: 73
[GPBoost] [Info] [GPBoost with gamma likelihood]: initscore=4.805531
[GPBoost] [Info] Start training from score 4.805531

19.7.7.7 Model Results

summary(gp_model_qb) # estimated random effects model

=====================================================
Covariance parameters (random effects):
                       Param.
Group_1                     0
Group_1_rand_coef_nb_1      0
Group_1_rand_coef_nb_2      0
-----------------------------------------------------
Additional parameters:
      Param.
shape 0.8186
=====================================================

gp_model_qb_importance <- gpboost::gpb.importance(gp_model_fit_qb)
gp_model_qb_importance

gpboost::gpb.plot.importance(gp_model_qb_importance)

Figure 19.8: Importance of Features (Predictors) in Tree Boosting Machine Learning Model.

19.7.7.8 Evaluate Accuracy of Model on Test Data

# Test Model on Test Data
pred_test_qb <- predict(
  gp_model_fit_qb,
  data = data_test_qb_matrix[,pred_vars_qb],
  group_data_pred = data_test_qb_matrix[,"gsis_id"],
  group_rand_coef_data_pred = cbind(
    data_test_qb_matrix[,"ageCentered20"],
    data_test_qb_matrix[,"ageCentered20Quadratic"]),
  predict_var = FALSE,
  pred_latent = FALSE)

y_pred_test_qb <- pred_test_qb[["response_mean"]] # if outcome is mean-centered, add mean(data_train_qb_matrix[,"fantasyPoints_lag"])

predictedVsActual <- data.frame(
  predictedPoints = y_pred_test_qb,
  actualPoints = data_test_qb_matrix[,"fantasyPoints_lag"]
)

predictedVsActual

petersenlab::accuracyOverall(
  predicted = predictedVsActual$predictedPoints,
  actual = predictedVsActual$actualPoints,
  dropUndefined = TRUE
)

19.7.7.9 Generate Predictions for Next Season

# Generate model predictions for next season
pred_nextYear_qb <- predict(
  gp_model_fit_qb,
  data = newData_qb_matrix[,pred_vars_qb],
  group_data_pred = newData_qb_matrix[,"gsis_id"],
  group_rand_coef_data_pred = cbind(
    newData_qb_matrix[,"ageCentered20"],
    newData_qb_matrix[,"ageCentered20Quadratic"]),
  predict_var = FALSE,
  pred_latent = FALSE)

newData_qb$fantasyPoints_lag <- pred_nextYear_qb$response_mean

# Merge with player names
newData_qb <- left_join(
  newData_qb,
  nfl_playerIDs %>% select(gsis_id, name),
  by = "gsis_id"
)

newData_qb %>% 
  arrange(-fantasyPoints_lag) %>% 
  select(name, fantasyPoints_lag, fantasyPoints)

19.8 Conclusion

19.9 Session Info

sessionInfo()

R version 4.5.1 (2025-06-13)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

time zone: UTC
tzcode source: system (glibc)

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] glmnet_4.1-9       Matrix_1.7-3       lubridate_1.9.4    forcats_1.0.0     
 [5] stringr_1.5.1      readr_2.1.5        tidyverse_2.0.0    gpboost_1.5.8     
 [9] R6_2.6.1           LongituRF_0.9      yardstick_1.3.2    workflowsets_1.1.1
[13] workflows_1.2.0    tune_1.3.0         tidyr_1.3.1        tibble_3.3.0      
[17] rsample_1.3.0      recipes_1.3.1      purrr_1.0.4        parsnip_1.3.2     
[21] modeldata_1.4.0    infer_1.0.9        ggplot2_3.5.2      dplyr_1.1.4       
[25] dials_1.4.0        scales_1.4.0       broom_1.0.8        tidymodels_1.3.0  
[29] powerjoin_0.1.0    missRanger_2.6.1   future_1.58.0      petersenlab_1.1.6 

loaded via a namespace (and not attached):
 [1] RColorBrewer_1.1-3   shape_1.4.6.1        rstudioapi_0.17.1   
 [4] jsonlite_2.0.0       magrittr_2.0.3       farver_2.1.2        
 [7] nloptr_2.2.1         rmarkdown_2.29       vctrs_0.6.5         
[10] minqa_1.2.8          base64enc_0.1-3      sparsevctrs_0.3.4   
[13] htmltools_0.5.8.1    Formula_1.2-5        parallelly_1.45.0   
[16] htmlwidgets_1.6.4    plyr_1.8.9           lifecycle_1.0.4     
[19] iterators_1.0.14     pkgconfig_2.0.3      fastmap_1.2.0       
[22] rbibutils_2.3        digest_0.6.37        colorspace_2.1-1    
[25] furrr_0.3.1          Hmisc_5.2-3          labeling_0.4.3      
[28] latex2exp_0.9.6      randomForest_4.7-1.2 RJSONIO_2.0.0       
[31] timechange_0.3.0     compiler_4.5.1       withr_3.0.2         
[34] htmlTable_2.4.3      backports_1.5.0      DBI_1.2.3           
[37] psych_2.5.6          MASS_7.3-65          lava_1.8.1          
[40] tools_4.5.1          pbivnorm_0.6.0       foreign_0.8-90      
[43] ranger_0.17.0        future.apply_1.20.0  nnet_7.3-20         
[46] doFuture_1.1.1       glue_1.8.0           quadprog_1.5-8      
[49] nlme_3.1-168         grid_4.5.1           checkmate_2.3.2     
[52] cluster_2.1.8.1      reshape2_1.4.4       generics_0.1.4      
[55] gtable_0.3.6         tzdb_0.5.0           class_7.3-23        
[58] data.table_1.17.6    hms_1.1.3            foreach_1.5.2       
[61] pillar_1.10.2        mitools_2.4          splines_4.5.1       
[64] lhs_1.2.0            lattice_0.22-7       survival_3.8-3      
[67] FNN_1.1.4.1          tidyselect_1.2.1     mix_1.0-13          
[70] knitr_1.50           reformulas_0.4.1     gridExtra_2.3       
[73] stats4_4.5.1         xfun_0.52            hardhat_1.4.1       
[76] timeDate_4041.110    stringi_1.8.7        DiceDesign_1.10     
[79] yaml_2.3.10          boot_1.3-31          evaluate_1.0.4      
[82] codetools_0.2-20     cli_3.6.5            rpart_4.1.24        
[85] xtable_1.8-4         Rdpack_2.6.4         lavaan_0.6-19       
[88] Rcpp_1.0.14          globals_0.18.0       gower_1.0.2         
[91] GPfit_1.0-9          lme4_1.1-37          listenv_0.9.1       
[94] viridisLite_0.4.2    mvtnorm_1.3-3        ipred_0.9-15        
[97] prodlim_2025.04.28   rlang_1.1.6          mnormt_2.1.1        

Bengtsson, H. (2025). ‘future‘: Unified parallel and distributed processing in ‘R‘ for everyone. https://doi.org/10.32614/CRAN.package.future

Capitaine, L. (2020). LongituRF: Random forests for longitudinal data. https://doi.org/10.32614/CRAN.package.LongituRF

Corston, R., & Colman, A. M. (2000). A crash course in SPSS for Windows. Wiley-Blackwell.

Fabri, A. (2022). powerjoin: Extensions of dplyr and fuzzyjoin join functions. https://doi.org/10.32614/CRAN.package.powerjoin

Hu, J., & Szymczak, S. (2023). A review on longitudinal data analysis with random forest. Briefings in Bioinformatics, 24(2). https://doi.org/10.1093/bib/bbad002

Mayer, M. (2024). missRanger: Fast imputation of missing values. https://doi.org/10.32614/CRAN.package.missRanger

Sigrist, F., Gyger, T., & Kuendig, P. (2025). gpboost: Combining tree-boosting with gaussian process and mixed effects models. https://doi.org/10.32614/CRAN.package.gpboost